三角波

Triangle Waves

正弦波只包含一个频率分量，因此其频谱只有一个峰值。更复杂的波形，例如图 1-2中的小提琴录音，会产生具有多个峰值的离散傅里叶变换 (DFT)。本节我们将探讨波形与其频谱之间的关系。

A sinusoid contains only one frequency component, so its spectrum has only one peak. More complicated waveforms, like the violin recording in Figure 1-2, yield DFTs with many peaks. In this section we investigate the relationship between waveforms and their spectrums.

我先从三角波形开始，它就像是正弦波的直线版本。图 2-1显示了一个频率为 200 Hz 的三角波形。

I’ll start with a triangle waveform, which is like a straight-line version of a sinusoid. Figure 2-1 shows a triangle waveform with frequency 200 Hz.

图 2-1. 200 Hz 三角波信号的片段。

要生成三角波，您可以从以下代码开始thinkdsp.TriangleSignal：

To generate a triangle wave, you can start with a thinkdsp.TriangleSignal:

class TriangleSignal(Sinusoid):
    
    def evaluate(self, ts):
        周期数 = self.freq * ts + self.offset / PI2
        frac, _ = np.modf(cycles)
        ys = np.abs(frac - 0.5)
        ys = normalize(unbias(ys), self.amp)
        返回 ys

class TriangleSignal(Sinusoid):
    
    def evaluate(self, ts):
        cycles = self.freq * ts + self.offset / PI2
        frac, _ = np.modf(cycles)
        ys = np.abs(frac - 0.5)
        ys = normalize(unbias(ys), self.amp)
        return ys

TriangleSignal它继承__init__自Sinusoid，因此它接受相同的参数：freq，amp和offset。

TriangleSignal inherits __init__ from Sinusoid, so it takes the same arguments: freq, amp, and offset.

唯一的区别在于evaluate……正如我们之前看到的，ts是我们想要评估信号的采样时间序列。

The only difference is evaluate. As we saw before, ts is the sequence of sample times where we want to evaluate the signal.

生成三角波的方法有很多种。细节并不重要，但基本evaluate原理如下：

There are many ways to generate a triangle wave. The details are not important, but here’s how evaluate works:

1. cycles是自开始时间以来的循环次数。将循环np.modf次数拆分为小数部分（存储在 中）和整数部分（忽略）^。1frac

2. cycles is the number of cycles since the start time. np.modf splits the number of cycles into the fraction part, stored in frac, and the integer part, which is ignored.¹

3. frac这是一个以给定频率从 0 渐变到 1 的序列。减去 0.5 得到介于 -0.5 和 0.5 之间的值。取绝对值得到一个在 0.5 和 0 之间呈锯齿状的波形。

4. frac is a sequence that ramps from 0 to 1 with the given frequency. Subtracting 0.5 yields values between –0.5 and 0.5. Taking the absolute value yields a waveform that zigzags between 0.5 and 0.

5. unbias将波形向下移动，使其中心位于 0；然后normalize将其缩放到给定的振幅amp。

6. unbias shifts the waveform down so it is centered at 0; then normalize scales it to the given amplitude, amp.

以下是生成图 2-1 的代码：

Here’s the code that generates Figure 2-1:

signal = thinkdsp.TriangleSignal(200)
signal.plot()

signal = thinkdsp.TriangleSignal(200)
signal.plot()

接下来，我们可以用它Signal来制作一个Wave，并用它Wave来制作一个Spectrum：

Next we can use the Signal to make a Wave, and use the Wave to make a Spectrum:

wave = signal.make_wave(duration=0.5, framerate=10000)
spectrum = wave.make_spectrum()
频谱图()

wave = signal.make_wave(duration=0.5, framerate=10000)
spectrum = wave.make_spectrum()
spectrum.plot()

图 2-2显示了结果的两个视图；右侧视图经过缩放，以便更清晰地显示谐波。正如预期的那样，最高峰值位于基频 200 Hz 处，并且在 200 的整数倍谐波频率处也存在其他峰值。

Figure 2-2 shows two views of the result; the view on the right is scaled to show the harmonics more clearly. As expected, the highest peak is at the fundamental frequency, 200 Hz, and there are additional peaks at harmonic frequencies, which are integer multiples of 200.

图 2-2. 200 Hz 三角波信号的频谱，以两个纵坐标刻度显示。右侧版本截去了基频，以便更清晰地显示谐波。

但令人惊讶的是，偶数倍频处没有峰值：400、800 等。三角波的谐波都是基频的奇数倍，在本例中为 600、1000、1400 等。

But one surprise is that there are no peaks at the even multiples: 400, 800, etc. The harmonics of a triangle wave are all odd multiples of the fundamental frequency, in this example 600, 1000, 1400, etc.

该频谱的另一个特征是谐波的振幅与频率之间的关系。谐波的振幅与频率的平方成正比下降。例如，前两个谐波（200 Hz 和 600 Hz）的频率比为 3，振幅比约为 9。接下来的两个谐波（600 Hz 和 1000 Hz）的频率比为 1.7，振幅比约为 1.7 。这种关系被称为谐波结构。

Another feature of this spectrum is the relationship between the amplitude and frequency of the harmonics. Their amplitude drops off in proportion to frequency squared. For example, the frequency ratio of the first two harmonics (200 and 600 Hz) is 3, and the amplitude ratio is approximately 9. The frequency ratio of the next two harmonics (600 and 1000 Hz) is 1.7, and the amplitude ratio is approximately . This relationship is called the harmonic structure.

方波

Square Waves

thinkdsp它还提供了SquareSignal一个表示方波信号的函数。以下是该类的定义：

thinkdsp also provides SquareSignal, which represents a square signal. Here’s the class definition:

class SquareSignal(Sinusoid):
    
    def evaluate(self, ts):
        周期数 = self.freq * ts + self.offset / PI2
        frac, _ = np.modf(cycles)
        ys = self.amp * np.sign(unbias(frac))
        返回 ys

class SquareSignal(Sinusoid):
    
    def evaluate(self, ts):
        cycles = self.freq * ts + self.offset / PI2
        frac, _ = np.modf(cycles)
        ys = self.amp * np.sign(unbias(frac))
        return ys

就像TriangleSignal，SquareSignal继承__init__自Sinusoid，所以它接受相同的参数。

Like TriangleSignal, SquareSignal inherits __init__ from Sinusoid, so it takes the same parameters.

方法evaluate类似。同样，cycles是自开始时间以来的周期数，frac是小数部分，每个周期从 0 递增到 1。

And the evaluate method is similar. Again, cycles is the number of cycles since the start time, and frac is the fractional part, which ramps from 0 to 1 each period.

unbias进行偏移frac，使其从 -0.5 逐渐过渡到 0.5，然后np.sign将负值映射到 -1，正值映射到 1。乘以得到一个在和amp之间跳跃的方波。-ampamp

unbias shifts frac so it ramps from –0.5 to 0.5, then np.sign maps the negative values to –1 and the positive values to 1. Multiplying by amp yields a square wave that jumps between -amp and amp.

图 2-3显示了频率为 100 Hz 的方波的三个周期，图 2-4显示了它的频谱。

Figure 2-3 shows three periods of a square wave with frequency 100 Hz, and Figure 2-4 shows its spectrum.

与三角波类似，方波也只包含奇次谐波，因此在 300、500 和 700 Hz 等频率处会出现峰值。但方波谐波的振幅衰减速度较慢。具体来说，振幅的衰减与频率成正比（而非与频率的平方成正比）。

Like a triangle wave, the square wave contains only odd harmonics, which is why there are peaks at 300, 500, and 700 Hz, etc. But the amplitude of the harmonics drops off more slowly. Specifically, amplitude drops in proportion to frequency (not frequency squared).

本章末尾的练习将使您有机会探索其他波形和其他谐波结构。

The exercises at the end of this chapter give you a chance to explore other waveforms and other harmonic structures.

图 2-3. 100 Hz 方波信号的片段。

图 2-4. 100 Hz 方波信号的频谱。

别名

Aliasing

我得坦白一件事。上一节的例子都是我精心挑选的，就是为了避免让你们感到困惑。但现在，是时候让你们感到困惑了。

I have a confession. I chose the examples in the previous section carefully to avoid showing you something confusing. But now it’s time to get confused.

图 2-5显示了频率为 1100 Hz、采样率为每秒 10,000 帧的三角波频谱。同样，右侧视图已按比例缩放以显示谐波。

Figure 2-5 shows the spectrum of a triangle wave at 1100 Hz, sampled at 10,000 frames per second. Again, the view on the right is scaled to show the harmonics.

图 2-5.采样频率为 1100 Hz、采样率为 10,000 帧/秒的三角波信号的频谱。右侧视图按比例缩放以显示谐波。

该波的谐波频率应为 3300、5500、7700 和 9900 Hz。图中，1100 和 3300 Hz 处有峰值，符合预期，但第三个峰值在 4500 Hz，而不是 5500 Hz。第四个峰值在 2300 Hz，而不是 7700 Hz。仔细观察，本应在 9900 Hz 的峰值实际上在 100 Hz。这是怎么回事？

The harmonics of this wave should be at 3300, 5500, 7700, and 9900 Hz. In the figure, there are peaks at 1100 and 3300 Hz, as expected, but the third peak is at 4500, not 5500 Hz. The fourth peak is at 2300, not 7700 Hz. And if you look closely, the peak that should be at 9900 is actually at 100 Hz. What’s going on?

问题在于，当你在离散的时间点评估信号时，你会丢失样本间发生的信息。对于低频分量来说，这不是问题，因为每个周期内有很多样本。

The problem is that when you evaluate the signal at discrete points in time, you lose information about what happened between samples. For low-frequency components, that’s not a problem, because you have lots of samples per period.

但如果你以 5000 Hz 的频率和每秒 10000 帧的速度对信号进行采样，那么每个周期你只有两个样本。这勉强够用，但如果频率更高，那就不够了。

But if you sample a signal at 5000 Hz with 10,000 frames per second, you only have two samples per period. That turns out to be enough, just barely, but if the frequency is higher, it’s not.

为了了解原因，我们生成频率分别为 4500 Hz 和 5500 Hz 的余弦信号，并以每秒 10,000 帧的速度对它们进行采样：

To see why, let’s generate cosine signals at 4500 and 5500 Hz, and sample them at 10,000 frames per second:

帧率 = 10000

signal = thinkdsp.CosSignal(4500)
持续时间 = 信号周期 * 5
segment = signal.make_wave(duration, framerate=framerate)
segment.plot()

signal = thinkdsp.CosSignal(5500)
segment = signal.make_wave(duration, framerate=framerate)
segment.plot()

framerate = 10000

signal = thinkdsp.CosSignal(4500)
duration = signal.period*5
segment = signal.make_wave(duration, framerate=framerate)
segment.plot()

signal = thinkdsp.CosSignal(5500)
segment = signal.make_wave(duration, framerate=framerate)
segment.plot()

图 2-6显示了结果。我Signal用细灰线绘制了 s 值，用垂直线绘制了样本值，以便于比较两个Waves 值。问题应该很明显：尽管Signals 值不同，但Waves 值却相同！

Figure 2-6 shows the result. I plotted the Signals with thin gray lines and the samples using vertical lines, to make it easier to compare the two Waves. The problem should be clear: even though the Signals are different, the Waves are identical!

图 2-6.采样频率分别为 4500 Hz 和 5500 Hz 的余弦信号，采样率为每秒 10,000 帧。信号本身不同，但采样点相同。

当我们以每秒 10,000 帧的速率对 5500 Hz 的信号进行采样时，结果与 4500 Hz 的信号无法区分。同样，7700 Hz 的信号与 2300 Hz 的信号无法区分，9900 Hz 的信号与 100 Hz 的信号无法区分。

When we sample a 5500 Hz signal at 10,000 frames per second, the result is indistinguishable from a 4500 Hz signal. For the same reason, a 7700 Hz signal is indistinguishable from 2300 Hz, and a 9900 Hz signal is indistinguishable from 100 Hz.

这种效应称为混叠，因为当对高频信号进行采样时，它看起来像是一个低频信号。

This effect is called aliasing because when the high-frequency signal is sampled, it appears to be a low-frequency signal.

在这个例子中，我们能测量的最高频率是 5000 Hz，也就是采样率的一半。高于 5000 Hz 的频率会被折叠回 5000 Hz 以下，这就是为什么这个阈值有时被称为“折叠频率”。它有时也被称为奈奎斯特频率。参见http://en.wikipedia.org/wiki/Nyquist_frequency。

In this example, the highest frequency we can measure is 5000 Hz, which is half the sampling rate. Frequencies above 5000 Hz are folded back below 5000 Hz, which is why this threshold is sometimes called the “folding frequency”. It is sometimes also called the Nyquist frequency. See http://en.wikipedia.org/wiki/Nyquist_frequency.

如果混叠频率低于零，折叠模式会继续。例如，1100 Hz 三角波的五次谐波频率为 12100 Hz。它在 5000 Hz 处折叠后，会出现在 -2100 Hz 处，但它会在 0 Hz 处再次折叠，回到 2100 Hz。实际上，您可以在图 2-4中看到 2100 Hz 处的小峰，以及 4300 Hz 处的下一个峰。

The folding pattern continues if the aliased frequency goes below zero. For example, the fifth harmonic of the 1100 Hz triangle wave is at 12,100 Hz. Folded at 5000 Hz, it would appear at –2100 Hz, but it gets folded again at 0 Hz, back to 2100 Hz. In fact, you can see a small peak at 2100 Hz in Figure 2-4, and the next one at 4300 Hz.

计算频谱

Computing the Spectrum

我们已经多次见过这种Wave方法make_spectrum。以下是它的实现方式（省略了一些细节，稍后会详细介绍）：

We have seen the Wave method make_spectrum several times. Here is the implementation (leaving out some details we’ll get to later):

from np.fft import rfft, rfftfreq

# Wave 类：
    def make_spectrum(self):
        n = len(self.ys)
        d = 1 / 自身帧率

        hs = rfft(self.ys)
        fs = rfftfreq(n, d)

        返回 Spectrum(hs, fs, self.framerate)

from np.fft import rfft, rfftfreq

# class Wave:
    def make_spectrum(self):
        n = len(self.ys)
        d = 1 / self.framerate

        hs = rfft(self.ys)
        fs = rfftfreq(n, d)

        return Spectrum(hs, fs, self.framerate)

该参数self是一个Wave对象。n是波形中的样本数，d是帧速率的倒数，即样本之间的时间间隔。

The parameter self is a Wave object. n is the number of samples in the wave, and d is the inverse of the frame rate, which is the time between samples.

np.fft是 NumPy 模块，提供与快速傅里叶变换 (FFT) 相关的函数，FFT 是一种计算离散傅里叶变换 (DFT) 的高效算法。

np.fft is the NumPy module that provides functions related to the Fast Fourier Transform (FFT), which is an efficient algorithm that computes the Discrete Fourier Transform (DFT).

make_spectrum这里使用rfft`real FFT`，因为它Wave包含的是实数，而不是复数。稍后我们将看到完整的 FFT，它可以处理复信号（参见“实信号的 DFT”）。`real FFT` 的结果rfft，我称之为hs`ff_a`，是一个 NumPy 复数数组，表示波中每个频率分量的振幅和相位偏移。

make_spectrum uses rfft, which stands for “real FFT”, because the Wave contains real values, not complex. Later we’ll see the full FFT, which can handle complex signals (see “DFT of Real Signals”). The result of rfft, which I call hs, is a NumPy array of complex numbers that represents the amplitude and phase offset of each frequency component in the wave.

结果rfftfreq（我称之为fs）是一个数组，其中包含与对应的频率hs。

The result of rfftfreq, which I call fs, is an array that contains frequencies corresponding to the hs.

为了理解复数中的值hs，可以考虑以下两种思考复数的方式：

To understand the values in hs, consider these two ways to think about complex numbers:

• 复数是实部和虚部之和，通常记为i = (i - i) ，其中i是虚数单位。你可以把x和y看作是笛卡尔坐标。

• A complex number is the sum of a real part and an imaginary part, often written , where i is the imaginary unit . You can think of x and y as Cartesian coordinates.

• 复数也可以表示为模和复指数的乘积，即 A = ϕ ，其中A是模，ϕ是以弧度为单位的角度，也称为“辐角”。你可以把A和ϕ看作是极坐标。

• A complex number is also the product of a magnitude and a complex exponential, , where A is the magnitude and ϕ is the angle in radians, also called the “argument”. You can think of A and ϕ as polar coordinates.

每个值hs对应一个频率分量：其幅度与对应分量的振幅成正比；其角度为相位偏移。

Each value in hs corresponds to a frequency component: its magnitude is proportional to the amplitude of the corresponding component; its angle is the phase offset.

该类Spectrum提供两个只读属性 `m`amps和 ` anglesθ`，它们返回表示物体的幅值和角度的 NumPy 数组hs。绘制Spectrum物体时，我们通常绘制 `m`与amps`θ` 的关系图。有时绘制 `θ`与 `θ` 的关系fs图也很有用。anglesfs

The Spectrum class provides two read-only properties, amps and angles, which return NumPy arrays representing the magnitudes and angles of the hs. When we plot a Spectrum object, we usually plot amps versus fs. Sometimes it is also useful to plot angles versus fs.

尽管查看实部和虚部可能很诱人hs，但你几乎永远不需要这样做。我建议你将DFT视为一个振幅和相位偏移的向量，它们恰好以复数的形式编码。

Although it might be tempting to look at the real and imaginary parts of hs, you will almost never need to. I encourage you to think of the DFT as a vector of amplitudes and phase offsets that happen to be encoded in the form of complex numbers.

要修改 a Spectrum，您可以hs直接访问它。例如：

To modify a Spectrum, you can access the hs directly. For example:

spectrum.hs *= 2
spectrum.hs[spectrum.fs > cutoff] = 0

spectrum.hs *= 2
spectrum.hs[spectrum.fs > cutoff] = 0

第一行将矩阵的元素乘以 2，这使得所有分量的振幅都加倍。第二行仅将频率超过某个截止频率hs的元素设为 0 。hs

The first line multiplies the elements of hs by 2, which doubles the amplitudes of all components. The second line sets to 0 only the elements of hs where the corresponding frequency exceeds some cutoff frequency.

但Spectrum也提供了执行这些操作的方法：

But Spectrum also provides methods to perform these operations:

spectrum.scale(2)
spectrum.low_pass(截止频率)

spectrum.scale(2)
spectrum.low_pass(cutoff)

您可以访问http://greenteapress.com/thinkdsp.html阅读这些方法及其他方法的文档。

You can read the documentation of these methods and others at http://greenteapress.com/thinkdsp.html.

到目前为止，你应该对Signal` std::string` Wave、`std::string` 和 ` Spectrumstd::string` 类的工作原理有了更清晰的了解，但我还没有解释快速傅里叶变换的工作原理。这还需要几章的时间来讲解。

At this point you should have a better idea of how the Signal, Wave, and Spectrum classes work, but I have not explained how the Fast Fourier Transform works. That will take a few more chapters.

练习

Exercises

这些练习的答案在chap02soln.ipynb……

Solutions to these exercises are in chap02soln.ipynb.

线性啁啾

Linear Chirp

我们将从啁啾信号开始，这是一种频率可变的信号。thinkdsp它提供了一个称为正弦波的Signal信号Chirp，该信号可以线性地扫描一系列频率。

We’ll start with a chirp, which is a signal with variable frequency. thinkdsp provides a Signal called Chirp that makes a sinusoid that sweeps linearly through a range of frequencies.

这里有一个例子，它的频率范围从 220 赫兹扫到 880 赫兹，也就是从 A3 到 A5 的两个八度：

Here’s an example that sweeps from 220 to 880 Hz, which is two octaves from A3 to A5:

signal = thinkdsp.Chirp(start=220, end=880)
wave = signal.make_wave()

signal = thinkdsp.Chirp(start=220, end=880)
wave = signal.make_wave()

图 3-1显示了该波在起始、中间和结束附近的部分波形。显然，频率正在增加。

Figure 3-1 shows segments of this wave near the beginning, middle, and end. It’s clear that the frequency is increasing.

图 3-1.啁啾波形在开始、中间和结束附近。

在继续之前，我们先来看看它Chirp是如何实现的。以下是类定义：

Before we go on, let’s see how Chirp is implemented. Here is the class definition:

class Chirp(Signal):
    
    def __init__(self, start=440, end=880, amp=1.0):
        self.start = start
        self.end = end
        self.amp = amp

class Chirp(Signal):
    
    def __init__(self, start=440, end=880, amp=1.0):
        self.start = start
        self.end = end
        self.amp = amp

start和end分别是啁啾声开始和结束时的频率（单位：赫兹）。amp是振幅。

start and end are the frequencies, in Hz, at the start and end of the chirp. amp is amplitude.

以下是评估信号的函数：

Here is the function that evaluates the signal:

def evaluate(self, ts):
    freqs = np.linspace(self.start, self.end, len(ts)-1)
    返回 self._evaluate(ts, freqs)

def evaluate(self, ts):
    freqs = np.linspace(self.start, self.end, len(ts)-1)
    return self._evaluate(ts, freqs)

ts是信号需要评估的时间点序列；为了简化这个函数，我假设它们是等间隔的。

ts is the sequence of points in time where the signal should be evaluated; to keep this function simple, I assume they are equally spaced.

如果 的长度为tsn ，你可以把它看作是一系列时间间隔。为了计算每个时间间隔内的频率，我使用，它返回一个介于和 之间的np.linspaceNumPy 值数组。startend

If the length of ts is n, you can think of it as a sequence of intervals of time. To compute the frequency during each interval, I use np.linspace, which returns a NumPy array of values between start and end.

_evaluate是一个私有方法，用于完成其余的数学运算¹ :

_evaluate is a private method that does the rest of the math¹:

def _evaluate(self, ts, freqs):
    dts = np.diff(ts)
    dphis = PI2 * freqs * dts
    phases = np.cumsum(dphis)
    phases = np.insert(phases, 0, 0)
    ys = self.amp * np.cos(phases)
    返回 ys

def _evaluate(self, ts, freqs):
    dts = np.diff(ts)
    dphis = PI2 * freqs * dts
    phases = np.cumsum(dphis)
    phases = np.insert(phases, 0, 0)
    ys = self.amp * np.cos(phases)
    return ys

np.diff计算集合中相邻元素之间的差值ts，并返回每个间隔的长度（以秒为单位）。如果集合中的元素ts等距分布，dts则所有值都相同。

np.diff computes the difference between adjacent elements of ts, returning the length of each interval in seconds. If the elements of ts are equally spaced, the dts values are all the same.

下一步是确定每个时间间隔内相位的变化量。在“信号对象”一节中，我们看到，当频率恒定时，相位ϕ随时间线性增加：

The next step is to figure out how much the phase changes during each interval. In “Signal Objects” we saw that when frequency is constant, the phase, ϕ, increases linearly over time:

当频率是时间的函数时，在短时间间隔Δt内的相位变化为：

When frequency is a function of time, the change in phase during a short time interval, Δt, is:

在 Python 中，由于freqs`x`包含dts时间间隔，我们可以这样写：

In Python, since freqs contains and dts contains the time intervals, we can write:

dphis = PI2 * freqs * dts

dphis = PI2 * freqs * dts

现在，由于dphis包含了相位变化，我们可以通过将每个时间步长的变化相加来得到总相位：

Now, since dphis contains the changes in phase, we can get the total phase at each timestep by adding up the changes:

phases = np.cumsum(dphis)
phases = np.insert(phases, 0, 0)

phases = np.cumsum(dphis)
phases = np.insert(phases, 0, 0)

np.cumsum计算累积和，这几乎是我们所想要的，但它不是从 0 开始的。所以我会np.insert在开头加上一个 0。

np.cumsum computes the cumulative sum, which is almost what we want, but it doesn’t start at 0. So I use np.insert to add a 0 at the beginning.

结果是一个 NumPy 数组，其中i第 n 个元素包含第 n 个区间末尾的前i几项之和dphis；也就是说，是第 n 个区间末尾的总相位i。最后，np.cos计算波的振幅作为相位的函数（记住，相位以弧度表示）。

The result is a NumPy array where the ith element contains the sum of the first i terms from dphis; that is, the total phase at the end of the ith interval. Finally, np.cos computes the amplitude of the wave as a function of phase (remember that phase is expressed in radians).

如果你懂微积分，你可能会注意到当Δt趋于无穷大时，极限为：

If you know calculus, you might notice that the limit as Δt gets small is:

两边同时除以dt得到：

Dividing through by dt yields:

换句话说，频率是相位的导数。反之，相位是频率的积分。过去我们试图cumsum从频率推导出相位时，实际上是在近似积分。

In other words, frequency is the derivative of phase. Conversely, phase is the integral of frequency. When we used cumsum to go from frequency to phase, we were approximating integration.

指数啁啾

Exponential Chirp

当你聆听这段啁啾声时，你可能会注意到它的音调起初上升得很快，然后逐渐减慢。这段啁啾声跨越了两个八度，但跨越第一个八度只需要2/3秒，而跨越第二个八度则需要两倍的时间。

When you listen to this chirp, you might notice that the pitch rises quickly at first and then slows down. The chirp spans two octaves, but it only takes 2/3 s to span the first octave, and twice as long to span the second.

原因在于我们对音高的感知取决于频率的对数。因此，我们听到的两个音符之间的音程取决于它们频率的比值，而不是频率差。“音程”是音乐术语，指的是两个音高之间感知到的差异。

The reason is that our perception of pitch depends on the logarithm of frequency. As a result, the interval we hear between two notes depends on the ratio of their frequencies, not the difference. “Interval” is the musical term for the perceived difference between two pitches.

例如，八度音程是指两个音高之比为2的音程。因此，220赫兹到440赫兹的音程是一个八度音程，440赫兹到880赫兹的音程也是一个八度音程。频率差异更大，但比例相同。

For example, an octave is an interval where the ratio of two pitches is 2. So the interval from 220 to 440 Hz is one octave and the interval from 440 to 880 Hz is also one octave. The difference in frequency is bigger, but the ratio is the same.

因此，如果频率呈线性增加，就像线性啁啾声一样，则感知到的音调呈对数增加。

As a result, if frequency increases linearly, as in a linear chirp, the perceived pitch increases logarithmically.

如果希望感知到的音调线性增加，频率就必须呈指数增长。具有这种形状的信号称为指数啁啾信号。

If you want the perceived pitch to increase linearly, the frequency has to increase exponentially. A signal with that shape is called an exponential chirp.

以下是ExpoChirp类定义：

Here’s the ExpoChirp class definition:

class ExpoChirp(Chirp):
    
    def evaluate(self, ts):
        start, end = np.log10(self.start), np.log10(self.end)
        freqs = np.logspace(start, end, len(ts)-1)
        返回 self._evaluate(ts, freqs)

class ExpoChirp(Chirp):
    
    def evaluate(self, ts):
        start, end = np.log10(self.start), np.log10(self.end)
        freqs = np.logspace(start, end, len(ts)-1)
        return self._evaluate(ts, freqs)

np.linspace此版本evaluate使用，而不是np.logspace，它创建了一系列对数等距分布的频率，这意味着它们呈指数增长。

Instead of np.linspace, this version of evaluate uses np.logspace, which creates a series of frequencies whose logarithms are equally spaced, which means that they increase exponentially.

就是这样；其他一切都和之前一样Chirp。以下是生成该模型的代码：

That’s it; everything else is the same as Chirp. Here’s the code that makes one:

signal = thinkdsp.ExpoChirp(start=220, end=880)
wave = signal.make_wave(duration=1)

signal = thinkdsp.ExpoChirp(start=220, end=880)
wave = signal.make_wave(duration=1)

您可以收听示例chap03.ipynb并比较线性啁啾和指数啁啾。

You can listen to the examples in chap03.ipynb and compare the linear and exponential chirps.

啁啾声的频谱

Spectrum of a Chirp

你觉得计算一个啁啾信号的频谱会发生什么？这里有一个例子，它构建了一个持续一秒、一个八度的啁啾信号及其频谱：

What do you think happens if you compute the spectrum of a chirp? Here’s an example that constructs a one-second, one-octave chirp and its spectrum:

signal = thinkdsp.Chirp(start=220, end=440)
wave = signal.make_wave(duration=1)
spectrum = wave.make_spectrum()

signal = thinkdsp.Chirp(start=220, end=440)
wave = signal.make_wave(duration=1)
spectrum = wave.make_spectrum()

图 3-2显示了结果。频谱在 220 至 440 Hz 的每个频率上都有成分，其变化看起来有点像索伦之眼（参见http://en.wikipedia.org/wiki/Sauron）。

Figure 3-2 shows the result. The spectrum has components at every frequency from 220 to 440 Hz, with variations that look a little like the Eye of Sauron (see http://en.wikipedia.org/wiki/Sauron).

图 3-2.一秒、一个八度音程的啁啾信号的频谱。

在 220 至 440 Hz 之间，频谱大致平坦，这表明信号在该范围内每个频率停留的时间相等。基于此，你应该能够推测出指数啁啾信号的频谱是什么样的。

The spectrum is approximately flat between 220 and 440 Hz, which indicates that the signal spends equal time at each frequency in this range. Based on that observation, you should be able to guess what the spectrum of an exponential chirp looks like.

频谱图可以提供信号结构的一些线索，但它掩盖了频率和时间之间的关系。例如，我们无法通过观察这个频谱图来判断频率是升高了还是降低了，或者两者兼有。

The spectrum gives hints about the structure of the signal, but it obscures the relationship between frequency and time. For example, we cannot tell by looking at this spectrum whether the frequency went up or down, or both.

频谱图

Spectrogram

为了恢复频率和时间之间的关系，我们可以将啁啾信号分成若干段，并绘制每一段的频谱。这种结果称为短时傅里叶变换（STFT）。

To recover the relationship between frequency and time, we can break the chirp into segments and plot the spectrum of each segment. The result is called a Short-Time Fourier Transform (STFT).

短时傅里叶变换 (STFT) 的可视化方法有很多种，但最常见的是频谱图，它以时间轴为横坐标，频率轴为纵坐标。频谱图中的每一列都显示了一小段信号的频谱，并用颜色或灰度来表示振幅。

There are several ways to visualize an STFT, but the most common is a spectrogram, which shows time on the x-axis and frequency on the y-axis. Each column in the spectrogram shows the spectrum of a short segment, using color or grayscale to represent amplitude.

例如，我将计算这个啁啾信号的频谱图：

As an example, I’ll compute the spectrogram of this chirp:

signal = thinkdsp.Chirp(start=220, end=440)
wave = signal.make_wave(duration=1, framerate=11025)

signal = thinkdsp.Chirp(start=220, end=440)
wave = signal.make_wave(duration=1, framerate=11025)

Wave提供make_spectrogram，它返回一个Spectrogram对象：

Wave provides make_spectrogram, which returns a Spectrogram object:

频谱图 = wave.make_spectrogram(seg_length=512)
光谱图.plot(high=700)

spectrogram = wave.make_spectrogram(seg_length=512)
spectrogram.plot(high=700)

seg_length是每个片段中的样本数。我选择 512 是因为当样本数为 2 的幂时，FFT 的效率最高。

seg_length is the number of samples in each segment. I chose 512 because the FFT is most efficient when the number of samples is a power of 2.

图 3-3显示了结果。x 轴表示时间，范围从 0 到 1 秒。y 轴表示频率，范围从 0 到 700 Hz。我截掉了频谱图的顶部部分；完整的范围达到 5512.5 Hz，这是帧速率的一半。

Figure 3-3 shows the result. The x-axis shows time from 0 to 1 seconds. The y-axis shows frequency from 0 to 700 Hz. I cut off the top part of the spectrogram; the full range goes to 5512.5 Hz, which is half of the frame rate.

图 3-3.一秒、一个八度音程的啁啾声的频谱图。

频谱图清晰地显示频率随时间线性增加。然而，请注意，每列中的峰值在2-3个单元格内呈现模糊现象。这种模糊反映了频谱图分辨率的局限性。

The spectrogram shows clearly that frequency increases linearly over time. However, notice that the peak in each column is blurred across 2–3 cells. This blurring reflects the limited resolution of the spectrogram.

加博尔极限

The Gabor Limit

频谱图的时间分辨率就是各个片段的持续时间，这对应于频谱图中单元格的宽度。由于每个片段包含 512 帧，而每秒有 11,025 帧，因此每个片段的持续时间约为 0.046 秒。

The time resolution of the spectrogram is the duration of the segments, which corresponds to the width of the cells in the spectrogram. Since each segment is 512 frames, and there are 11,025 frames per second, the duration of each segment is about 0.046 seconds.

频率分辨率是指频谱中各元素之间的频率范围，它对应于单元格的高度。512 帧数据可以得到 256 个频率分量，范围从 0 到 5512.5 Hz，因此分量之间的频率范围为 21.6 Hz。

The frequency resolution is the frequency range between elements in the spectrum, which corresponds to the height of the cells. With 512 frames, we get 256 frequency components over a range from 0 to 5512.5 Hz, so the range between components is 21.6 Hz.

更一般地，如果n是片段长度，则频谱包含个分量。如果帧速率为r，则频谱中的最大频率为。因此，时间分辨率为，频率分辨率为 。

More generally, if n is the segment length, the spectrum contains components. If the frame rate is r, the maximum frequency in the spectrum is . So the time resolution is and the frequency resolution is:

即。

which is .

理想情况下，我们希望时间分辨率小，以便观察频率的快速变化；我们也希望频率分辨率小，以便观察频率的微小变化。但两者不可兼得。请注意，时间分辨率是频率分辨率 的倒数。因此，如果其中一个变小，另一个就会变大。

Ideally we would like time resolution to be small, so we can see rapid changes in frequency. And we would like frequency resolution to be small so we can see small changes in frequency. But you can’t have both. Notice that time resolution, , is the inverse of frequency resolution, . So if one gets smaller, the other gets bigger.

例如，如果将片段长度加倍，频率分辨率会减半（这是好事），但时间分辨率会加倍（这是坏事）。即使提高帧速率也无济于事。虽然采样点数量增加了，但频率范围也随之增大。

For example, if you double the segment length, you cut frequency resolution in half (which is good), but you double time resolution (which is bad). Even increasing the frame rate doesn’t help. You get more samples, but the range of frequencies increases at the same time.

这种权衡被称为Gabor 极限，它是这种时频分析的一个根本限制。

This tradeoff is called the Gabor limit and it is a fundamental limitation of this kind of time–frequency analysis.

窗口

Windowing

我们可以通过平滑段的开头和结尾之间的不连续性来减少泄漏，而实现这一目标的一种方法就是加窗。

We can reduce leakage by smoothing out the discontinuity between the beginning and end of the segment, and one way to do that is windowing.

“窗口”函数旨在将非周期性线段转换为看似周期性的线段。图 3-5（上图）显示了一个线段，其端点与起点并不平滑连接。

A “window” is a function designed to transform a non-periodic segment into something that can pass for periodic. Figure 3-5 (top) shows a segment where the end does not connect smoothly to the beginning.

图 3-5（中间）展示了“汉明窗”，这是一种比较常见的窗函数。没有完美的窗函数，但有些窗函数在不同的应用场景下表现最佳，而汉明窗就是一种用途广泛的优秀窗函数。

Figure 3-5 (middle) shows a “Hamming window”, one of the more common window functions. No window function is perfect, but some can be shown to be optimal for different applications, and Hamming is a good, all-purpose window.

图 3-5（底部）显示了将窗口乘以原始信号的结果。当窗口值接近 1 时，信号保持不变；当窗口值接近 0 时，信号被衰减。由于窗口在两端逐渐变窄，因此信号段的末端与起始端平滑衔接。

Figure 3-5 (bottom) shows the result of multiplying the window by the original signal. Where the window is close to 1, the signal is unchanged. Where the window is close to 0, the signal is attenuated. Because the window tapers at both ends, the end of the segment connects smoothly to the beginning.

图 3-5.正弦曲线段（上）、汉明窗（中）以及该曲线段与汉明窗的乘积（下）。

以下是代码。Wave它提供window了应用汉明窗的功能：

Here’s what the code looks like. Wave provides window, which applies a Hamming window:

#class Wave：
    def window(self, window):
        self.ys *= window

#class Wave:
    def window(self, window):
        self.ys *= window

NumPy 提供了hamming一个函数，用于计算给定长度的汉明窗：

And NumPy provides hamming, which computes a Hamming window with a given length:

window = np.hamming(len(wave))
wave.window(window)

window = np.hamming(len(wave))
wave.window(window)

NumPy 提供了计算其他窗口函数的函数，包括bartlett、blackman、hanning和kaiser。本章末尾的练习之一要求你尝试使用这些其他窗口函数。

NumPy provides functions to compute other window functions, including bartlett, blackman, hanning, and kaiser. One of the exercises at the end of this chapter asks you to experiment with these other windows.

实施频谱图

Implementing Spectrograms

现在我们理解了窗口化，就可以理解其实现方式了make_spectrogram。以下是Wave计算频谱图的方法：

Now that we understand windowing, we can understand the implementation of make_spectrogram. Here is the Wave method that computes spectrograms:

#class Wave：
    def make_spectrogram(self, seg_length):
        window = np.hamming(seg_length)
        i, j = 0, seg_length
        步长 = 段长度 / 2

        spec_map = {}

        当 j < len(self.ys):
            segment = self.slice(i, j)
            segment.window(window)

            t = (segment.start + segment.end) / 2
            spec_map[t] = segment.make_spectrum()

            i += 步骤
            j += 步骤

        返回频谱图(spec_map, seg_length)

#class Wave:
    def make_spectrogram(self, seg_length):
        window = np.hamming(seg_length)
        i, j = 0, seg_length
        step = seg_length / 2

        spec_map = {}

        while j < len(self.ys):
            segment = self.slice(i, j)
            segment.window(window)

            t = (segment.start + segment.end) / 2
            spec_map[t] = segment.make_spectrum()

            i += step
            j += step

        return Spectrogram(spec_map, seg_length)

这是本书中最长的函数，所以如果你能处理好这个函数，你就能处理好任何函数。

This is the longest function in the book, so if you can handle this, you can handle anything.

参数self是一个Wave对象，seg_length表示每个片段中的样本数。

The parameter, self, is a Wave object. seg_length is the number of samples in each segment.

window是一个与分段长度相同的汉明窗。

window is a Hamming window with the same length as the segments.

i和j是用于从波形中选择片段的切片索引。step是片段之间的偏移量。由于step是 的一半seg_length，因此片段重叠一半。图 3-6显示了这些重叠窗口的样子。

i and j are the slice indices that select segments from the wave. step is the offset between segments. Since step is half of seg_length, the segments overlap by half. Figure 3-6 shows what these overlapping windows look like.

图 3-6.重叠的汉明窗。

spec_map是一个字典，它将时间戳映射到Spectrum。

spec_map is a dictionary that maps from a timestamp to a Spectrum.

在循环内部while，我们从波形中选择一个切片并应用窗口；然后我们构造一个Spectrum对象并将其添加到spec_map。每个片段的标称时间，t即中点。

Inside the while loop, we select a slice from the wave and apply the window; then we construct a Spectrum object and add it to spec_map. The nominal time of each segment, t, is the midpoint.

然后我们前进i，j并继续，只要j不超过末尾Wave。

Then we advance i and j, and continue as long as j doesn’t go past the end of the Wave.

最后，该方法构造并返回一个Spectrogram对象。以下是该类的定义：

Finally, the method constructs and returns a Spectrogram object. Here is the definition of the class:

class Spectrogram(object):
    def __init__(self, spec_map, seg_length):
        self.spec_map = spec_map
        self.seg_length = seg_length

class Spectrogram(object):
    def __init__(self, spec_map, seg_length):
        self.spec_map = spec_map
        self.seg_length = seg_length

与许多初始化方法一样，这个方法只是将参数存储为属性。

Like many init methods, this one just stores the parameters as attributes.

Spectrogram提供plot，生成以时间为 x 轴、频率为 y 轴的伪彩色图。

Spectrogram provides plot, which generates a pseudocolor plot with time along the x-axis and frequency along the y-axis.

这就是Spectrograms 的实现方式。

And that’s how Spectrograms are implemented.

练习

Exercises

这些练习的答案在chap03soln.ipynb……

Solutions to these exercises are in chap03soln.ipynb.

不相关噪声

Uncorrelated Noise

理解噪声最简单的方法就是生成噪声，而最容易生成的噪声类型是无相关均匀噪声（UU噪声）。“均匀”是指信号包含来自均匀分布的随机值；也就是说，该范围内的每个值出现的概率均等。“无相关”是指这些值是独立的；也就是说，知道其中一个值并不能提供关于其他值的任何信息。

The simplest way to understand noise is to generate it, and the simplest kind to generate is uncorrelated uniform noise (UU noise). “Uniform” means the signal contains random values from a uniform distribution; that is, every value in the range is equally likely. “Uncorrelated” means that the values are independent; that is, knowing one value provides no information about the others.

这是一个表示 UU 噪声的类：

Here’s a class that represents UU noise:

class UncorrelatedUniformNoise(_Noise):

    def evaluate(self, ts):
        ys = np.random.uniform(-self.amp, self.amp, len(ts))
        返回 ys

class UncorrelatedUniformNoise(_Noise):

    def evaluate(self, ts):
        ys = np.random.uniform(-self.amp, self.amp, len(ts))
        return ys

UncorrelatedUniformNoise继承自_Noise，而又继承自Signal。

UncorrelatedUniformNoise inherits from _Noise, which inherits from Signal.

与往常一样，`evaluate` 函数接受 `t`ts参数，即信号需要被评估的时间点。它使用np.random.uniform`u` 函数，该函数从均匀分布中生成值。在本例中，这些值介于 `0`-amp和 `1`之间amp。

As usual, the evaluate function takes ts, the times when the signal should be evaluated. It uses np.random.uniform, which generates values from a uniform distribution. In this example, the values are in the range between -amp and amp.

以下示例生成持续时间为 0.5 秒、采样率为每秒 11,025 次的 UU 噪声：

The following example generates UU noise with duration 0.5 seconds at 11,025 samples per second:

signal = thinkdsp.UncorrelatedUniformNoise()
wave = signal.make_wave(duration=0.5, framerate=11025)

signal = thinkdsp.UncorrelatedUniformNoise()
wave = signal.make_wave(duration=0.5, framerate=11025)

播放这段波形，听起来就像收音机换台时听到的静电噪音。图 4-1显示了波形的样子。不出所料，它看起来相当随机。

If you play this wave, it sounds like the static you hear if you tune a radio between channels. Figure 4-1 shows what the waveform looks like. As expected, it looks pretty random.

图 4-1.不相关均匀噪声的波形。

现在我们来看一下频谱：

Now let’s take a look at the spectrum:

spectrum = wave.make_spectrum()
spectrum.plot_power()

spectrum = wave.make_spectrum()
spectrum.plot_power()

Spectrum.plot_power它与 类似Spectrum.plot，区别在于它绘制的是功率而不是振幅。功率是振幅的平方。本章我将从振幅切换到功率，是因为在噪声的语境下，功率更为常用。

Spectrum.plot_power is similar to Spectrum.plot, except that it plots power instead of amplitude. Power is the square of amplitude. I am switching from amplitude to power in this chapter because it is more conventional in the context of noise.

图 4-2显示了结果。与信号一样，频谱看起来相当随机。事实上，它是随机的，但我们需要更精确地定义“随机”一词。关于噪声信号或其频谱，我们至少需要了解以下三点：

Figure 4-2 shows the result. Like the signal, the spectrum looks pretty random. In fact, it is random, but we have to be more precise about the word “random”. There are at least three things we might like to know about a noise signal or its spectrum:

分配

Distribution

随机信号的分布是指所有可能取值及其概率的集合。例如，在均匀噪声信号中，取值范围为 -1 到 1，且所有值的概率均相同。另一种情况是高斯噪声，其取值范围为负无穷到正无穷，但接近 0 的值出现概率最高，其概率随时间呈高斯分布（或称“钟形”曲线）递减。

The distribution of a random signal is the set of possible values and their probabilities. For example, in the uniform noise signal, the set of values is the range from –1 to 1, and all values have the same probability. An alternative is Gaussian noise, where the set of values is the range from negative to positive infinity, but values near 0 are the most likely, with probability that drops off according to the Gaussian or “bell” curve.

相关性

Correlation

信号中的每个值都与其他值独立吗？还是它们之间存在依赖关系？在 UU 噪声中，各个值是独立的。另一种情况是布朗噪声，其中每个值都是前一个值加上一个随机“阶跃”之和。因此，如果信号在某一时刻的值很高，我们预期它会保持高值；如果它很低，我们预期它会保持低值。

Is each value in the signal independent from the others, or are there dependencies between them? In UU noise, the values are independent. An alternative is Brownian noise, where each value is the sum of the previous value and a random “step”. So if the value of the signal is high at a particular point in time, we expect it to stay high, and if it is low, we expect it to stay low.

功率与频率的关系

Relationship between power and frequency

在UU噪声频谱中，所有频率的功率都来自同一分布；也就是说，所有频率的平均功率都相同。另一种选择是粉红噪声，其功率与频率成反比；也就是说，频率f处的功率来自均值与频率成正比的分布。

In the spectrum of UU noise, the power at all frequencies is drawn from the same distribution; that is, the average power is the same for all frequencies. An alternative is pink noise, where power is inversely related to frequency; that is, the power at frequency f is drawn from a distribution whose mean is proportional to .

图 4-2.不相关均匀噪声的功率谱。

综合频谱

Integrated Spectrum

对于 UU 噪声，我们可以通过观察积分频谱来更清楚地看到功率和频率之间的关系，积分频谱是频率f的函数，它显示了频谱中直至f 的累积功率。

For UU noise we can see the relationship between power and frequency more clearly by looking at the integrated spectrum, which is a function of frequency, f, that shows the cumulative power in the spectrum up to f.

Spectrum提供一种计算以下值的方法IntegratedSpectrum：

Spectrum provides a method that computes the IntegratedSpectrum:

def make_integrated_spectrum(self):
    cs = np.cumsum(self.power)
    cs /= cs[-1]
    返回 IntegratedSpectrum(cs, self.fs)

def make_integrated_spectrum(self):
    cs = np.cumsum(self.power)
    cs /= cs[-1]
    return IntegratedSpectrum(cs, self.fs)

self.power是一个包含每个频率功率的 NumPy 数组。np.cumsum计算功率的累积和。除以最后一个元素，将积分频谱归一化，使其范围从 0 到 1。

self.power is a NumPy array containing the power for each frequency. np.cumsum computes the cumulative sum of the powers. Dividing through by the last element normalizes the integrated spectrum so it runs from 0 to 1.

结果是一个类IntegratedSpectrum。以下是类定义：

The result is an IntegratedSpectrum. Here is the class definition:

class IntegratedSpectrum(object):    
    def __init__(self, cs, fs):
        self.cs = cs
        self.fs = fs

class IntegratedSpectrum(object):    
    def __init__(self, cs, fs):
        self.cs = cs
        self.fs = fs

例如Spectrum，IntegratedSpectrum提供了plot_power，因此我们可以计算并绘制积分光谱，如下所示：

Like Spectrum, IntegratedSpectrum provides plot_power, so we can compute and plot the integrated spectrum like this:

integ = spectrum.make_integrated_spectrum()
integ.plot_power()

integ = spectrum.make_integrated_spectrum()
integ.plot_power()

如图 4-3所示，结果是一条直线，表明所有频率的功率平均而言是恒定的。所有频率功率相等的噪声被称为白噪声，这是类比于光，因为所有可见频率的等量混合光是白色的。

The result, shown in Figure 4-3, is a straight line, which indicates that power at all frequencies is constant, on average. Noise with equal power at all frequencies is called white noise by analogy with light, because an equal mixture of light at all visible frequencies is white.

图 4-3.不相关均匀噪声的积分频谱。

布朗噪声

Brownian Noise

UU噪声是不相关的，这意味着每个值与其他值无关。另一种方法是布朗噪声，其中每个值都是前一个值加上一个随机“步长”的总和。

UU noise is uncorrelated, which means that each value does not depend on the others. An alternative is Brownian noise, in which each value is the sum of the previous value and a random “step”.

它之所以被称为“布朗运动”，是因为它与布朗运动类似。在布朗运动中，悬浮在流体中的粒子由于与流体之间看不见的相互作用而看似随机地运动。布朗运动通常用随机游走来描述，随机游走是一种数学模型，它描述了一条路径，其中步长之间的距离由随机分布表征。

It is called “Brownian” by analogy with Brownian motion, in which a particle suspended in a fluid moves apparently at random, due to unseen interactions with the fluid. Brownian motion is often described using a random walk, which is a mathematical model of a path where the distance between steps is characterized by a random distribution.

在一维随机游走中，粒子在每个时间步长内上下移动的距离都是随机的。粒子在任何时刻的位置都是之前所有时间步长位置的总和。

In a one-dimensional random walk, the particle moves up or down by a random amount at each timestep. The location of the particle at any point in time is the sum of all previous steps.

这一观察结果提示了一种生成布朗噪声的方法：生成不相关的随机步长，然后将它们累加起来。以下是实现该算法的类定义：

This observation suggests a way to generate Brownian noise: generate uncorrelated random steps and then add them up. Here is a class definition that implements this algorithm:

class BrownianNoise(_Noise):

    def evaluate(self, ts):
        dy = np.random.uniform(-1, 1, len(ts))
        ys = np.cumsum(dys)
        ys = normalize(unbias(ys), self.amp)
        返回 ys

class BrownianNoise(_Noise):

    def evaluate(self, ts):
        dys = np.random.uniform(-1, 1, len(ts))
        ys = np.cumsum(dys)
        ys = normalize(unbias(ys), self.amp)
        return ys

evaluate用于np.random.uniform生成不相关信号并np.cumsum计算其累积和。

evaluate uses np.random.uniform to generate an uncorrelated signal and np.cumsum to compute their cumulative sum.

由于总和很可能超出 -1 到 1 的范围，我们必须使用unbias将均值移至 0，以normalize获得所需的最大振幅。

Since the sum is likely to escape the range from –1 to 1, we have to use unbias to shift the mean to 0, and normalize to get the desired maximum amplitude.

以下代码用于生成BrownianNoise对象并绘制波形：

Here’s the code that generates a BrownianNoise object and plots the waveform:

signal = thinkdsp.BrownianNoise()
wave = signal.make_wave(duration=0.5, framerate=11025)
wave.plot()

signal = thinkdsp.BrownianNoise()
wave = signal.make_wave(duration=0.5, framerate=11025)
wave.plot()

图 4-4显示了结果。波形上下波动，但相邻数值之间存在明显的关联性。当振幅较高时，它往往保持较高水平，反之亦然。

Figure 4-4 shows the result. The waveform wanders up and down, but there is a clear correlation between successive values. When the amplitude is high, it tends to stay high, and vice versa.

图 4-4.布朗噪声波形。

如果将布朗噪声的频谱绘制成线性坐标图，如图4-5（左）所示，看起来并不复杂。几乎所有的能量都集中在最低频率上；高频分量则不可见。

If you plot the spectrum of Brownian noise on a linear scale, as in Figure 4-5 (left), it doesn’t look like much. Nearly all of the power is at the lowest frequencies; the higher frequency components are not visible.

图 4-5.布朗噪声的频谱，线性尺度（左）和双对数尺度（右）。

为了更清晰地观察频谱形状，我们可以用双对数坐标绘制功率和频率曲线。以下是代码：

To see the shape of the spectrum more clearly, we can plot power and frequency on a log-log scale. Here’s the code:

spectrum = wave.make_spectrum()
spectrum.plot_power(linewidth=1, alpha=0.5)
thinkplot.config(xscale='log', yscale='log')

spectrum = wave.make_spectrum()
spectrum.plot_power(linewidth=1, alpha=0.5)
thinkplot.config(xscale='log', yscale='log')

结果如图 4-5（右）所示。功率与频率之间的关系存在噪声，但大致呈线性关系。

The result is in Figure 4-5 (right). The relationship between power and frequency is noisy, but roughly linear.

Spectrum提供estimate_slope，它使用 SciPy 计算功率谱的最小二乘拟合：

Spectrum provides estimate_slope, which uses SciPy to compute a least squares fit to the power spectrum:

#光谱类

    def estimate_slope(self):
        x = np.log(self.fs[1:])
        y = np.log(self.power[1:])
        t = scipy.stats.linregress(x,y)
        返回 t

#class Spectrum

    def estimate_slope(self):
        x = np.log(self.fs[1:])
        y = np.log(self.power[1:])
        t = scipy.stats.linregress(x,y)
        return t

它舍弃了频谱的第一个分量，因为该分量对应于，并且未定义。

It discards the first component of the spectrum because this component corresponds to , and is undefined.

estimate_slope返回结果scipy.stats.linregress，该结果是一个包含估计斜率和截距、决定系数 ( R² ⁾、p 值和标准误差的对象。就我们的目的而言，我们只需要斜率。

estimate_slope returns the result from scipy.stats.linregress, which is an object that contains the estimated slope and intercept, coefficient of determination (R²), p-value, and standard error. For our purposes, we only need the slope.

对于布朗噪声，功率谱的斜率为 -2（我们将在第 9 章中看到原因），因此我们可以写出以下关系：

For Brownian noise, the slope of the power spectrum is –2 (we’ll see why in Chapter 9), so we can write this relationship:

其中P为功率，f为频率，k为直线的截距，但对于我们的目的而言，截距并不重要。两边同时取指数可得：

where P is power, f is frequency, and k is the intercept of the line, which is not important for our purposes. Exponentiating both sides yields:

其中K是e ^k，但这仍然不重要。更重要的是，功率与成正比，这是布朗噪声的特征。

where K is e^k, but still not important. More relevant is that power is proportional to , which is characteristic of Brownian noise.

布朗噪声也被称为红噪声，原因与白噪声被称为“白噪声”相同。如果将可见光与功率成正比的辐射混合，大部分功率将集中在频谱的低频端，也就是红色区域。布朗噪声有时也被称为“棕色噪声”，但我认为这容易造成混淆，所以我就不这么叫了。

Brownian noise is also called red noise, for the same reason that white noise is called “white”. If you combine visible light with power proportional to , most of the power will be at the low-frequency end of the spectrum, which is red. Brownian noise is also sometimes called “brown noise”, but I think that’s confusing, so I won’t use it.

粉红噪音

Pink Noise

对于红噪声，频率和功率之间的关系为：

For red noise, the relationship between frequency and power is:

指数 2 并没有什么特殊之处。更一般地，我们可以用任意指数β合成噪声：

There is nothing special about the exponent 2. More generally, we can synthesize noise with any exponent, β:

当功率在所有频率上都恒定时，结果为白噪声。当功率增大时，结果为红噪声。

When , power is constant at all frequencies, so the result is white noise. When the result is red noise.

当β介于0和2之间时，结果介于白噪声和红噪声之间，因此被称为粉红噪声。

When β is between 0 and 2, the result is between white and red noise, so it is called pink noise.

生成粉红噪声的方法有很多种。最简单的方法是先生成白噪声，然后应用所需指数的低通滤波器。以下thinkdsp提供了一个表示粉红噪声信号的类：

There are several ways to generate pink noise. The simplest is to generate white noise and then apply a low-pass filter with the desired exponent. thinkdsp provides a class that represents a pink noise signal:

class PinkNoise(_Noise):

    def __init__(self, amp=1.0, beta=1.0):
        self.amp = amp
        self.beta = beta

class PinkNoise(_Noise):

    def __init__(self, amp=1.0, beta=1.0):
        self.amp = amp
        self.beta = beta

amp是信号的期望幅度。beta是期望指数。PinkNoise提供make_wave，它生成Wave：

amp is the desired amplitude of the signal. beta is the desired exponent. PinkNoise provides make_wave, which generates a Wave:

def make_wave(self, duration=1, start=0, framerate=11025):
   信号 = UncorrelatedUniformNoise()
    wave = signal.make_wave(duration, start, framerate)
    spectrum = wave.make_spectrum()

    spectrum.pink_filter(beta=self.beta)

    wave2 = spectrum.make_wave()
    wave2.unbias()
    wave2.normalize(self.amp)
    返回波2

def make_wave(self, duration=1, start=0, framerate=11025):
   signal = UncorrelatedUniformNoise()
    wave = signal.make_wave(duration, start, framerate)
    spectrum = wave.make_spectrum()

    spectrum.pink_filter(beta=self.beta)

    wave2 = spectrum.make_wave()
    wave2.unbias()
    wave2.normalize(self.amp)
    return wave2

duration是波的长度，单位为秒。start是波的起始时间；包含起始时间是为了确保make_wave所有类型的信号都能使用相同的接口，但对于随机噪声，起始时间无关紧要。framerate是每秒采样数。

duration is the length of the wave in seconds. start is the start time of the wave; it is included so that make_wave has the same interface for all types of signal, but for random noise, start time is irrelevant. framerate is the number of samples per second.

make_wave生成一个白噪声波，计算其频谱，应用所需指数的滤波器，然后将滤波后的频谱转换回波形。最后，对波形进行去偏和归一化处理。

make_wave creates a white noise wave, computes its spectrum, applies a filter with the desired exponent, and then converts the filtered spectrum back to a wave. Then it unbiases and normalizes the wave.

Spectrum提供pink_filter：

Spectrum provides pink_filter:

def pink_filter(self, beta=1.0):
    denom = self.fs ** (beta/2.0)
    denom[0] = 1
    self.hs /= denom

def pink_filter(self, beta=1.0):
    denom = self.fs ** (beta/2.0)
    denom[0] = 1
    self.hs /= denom

pink_filter将频谱的每个元素除以。由于功率是幅度的平方，因此该操作将每个分量的功率除以f ^β。它将 处分量视为特例，部分原因是避免除以 0，部分原因是该元素代表信号的偏置，而我们无论如何都要将其设置为 0。

pink_filter divides each element of the spectrum by . Since power is the square of amplitude, this operation divides the power at each component by f^β. It treats the component at as a special case, partly to avoid dividing by 0 and partly because this element represents the bias of the signal, which we are going to set to 0 anyway.

图 4-6显示了所得波形。与布朗噪声类似，它的波形上下波动，似乎暗示着连续值之间存在相关性，但至少从视觉上看，它更偏向随机性。下一章我们将再次讨论这一观察结果，届时我将更精确地解释“相关性”和“更偏向随机性”的含义。

Figure 4-6 shows the resulting waveform. Like that of Brownian noise, it wanders up and down in a way that suggests correlation between successive values, but at least visually, it looks more random. In the next chapter we will come back to this observation and I will be more precise about what I mean by “correlation” and “more random”.

图 4-6.粉红噪声波形。

最后，图 4-7显示了白噪声、粉红噪声和红噪声在同一双对数坐标系下的频谱。图中可以明显看出指数β与频谱斜率之间的关系。

Finally, Figure 4-7 shows a spectrum for white, pink, and red noise on the same log-log scale. The relationship between the exponent, β, and the slope of the spectrum is apparent in this figure.

图 4-7.白噪声、粉红噪声和红噪声的频谱（对数-对数坐标）。

高斯噪声

Gaussian Noise

我们从不相关的均匀（UU）噪声开始，并证明，由于其频谱在所有频率上都具有相同的功率，因此平均而言，UU 噪声是白噪声。

We started with uncorrelated uniform (UU) noise and showed that, because its spectrum has equal power at all frequencies, on average, UU noise is white.

但人们谈论“白噪声”时，并不总是指UU噪声。事实上，他们更多时候指的是不相关的高斯（UG）噪声。

But when people talk about “white noise”, they don’t always mean UU noise. In fact, more often they mean uncorrelated Gaussian (UG) noise.

thinkdsp提供了一种 UG 噪声的实现方式：

thinkdsp provides an implementation of UG noise:

class UncorrelatedGaussianNoise(_Noise):

    def evaluate(self, ts):
        ys = np.random.normal(0, self.amp, len(ts))
        返回 ys

class UncorrelatedGaussianNoise(_Noise):

    def evaluate(self, ts):
        ys = np.random.normal(0, self.amp, len(ts))
        return ys

np.random.normal返回一个 NumPy 数组，其中包含服从高斯分布的值，在本例中均值为 0，标准差为 σ self.amp。理论上，取值范围从负无穷到正无穷，但我们预期约 99% 的值介于 -3 和 3 之间。

np.random.normal returns a NumPy array of values from a Gaussian distribution, in this case with mean 0 and standard deviation self.amp. In theory the range of values is from negative to positive infinity, but we expect about 99% of the values to be between –3 and 3.

UG噪声在很多方面与UU噪声相似。平均而言，其频谱在所有频率上的功率相等，因此UG噪声也是白噪声。它还有一个有趣的特性：UG噪声的频谱本身也是UG噪声。更准确地说，频谱的实部和虚部是不相关的高斯值。

UG noise is similar in many ways to UU noise. The spectrum has equal power at all frequencies, on average, so UG is also white. And it has one other interesting property: the spectrum of UG noise is also UG noise. More precisely, the real and imaginary parts of the spectrum are uncorrelated Gaussian values.

为了验证这一说法，我们可以生成 UG 噪声的频谱，然后生成“正态概率图”，这是一种用图形方式检验分布是否为高斯分布的方法：

To test that claim, we can generate the spectrum of UG noise and then generate a “normal probability plot”, which is a graphical way to test whether a distribution is Gaussian:

signal = thinkdsp.UncorrelatedGaussianNoise()
wave = signal.make_wave(duration=0.5, framerate=11025)
spectrum = wave.make_spectrum()

thinkstats2.NormalProbabilityPlot(spectrum.real)
thinkstats2.NormalProbabilityPlot(spectrum.imag)

signal = thinkdsp.UncorrelatedGaussianNoise()
wave = signal.make_wave(duration=0.5, framerate=11025)
spectrum = wave.make_spectrum()

thinkstats2.NormalProbabilityPlot(spectrum.real)
thinkstats2.NormalProbabilityPlot(spectrum.imag)

NormalProbabilityPlot本书的存储库中包含了该工具thinkstats2。如果您不熟悉正态概率图，可以在Think Stats的第 5 章中阅读相关内容，网址为http://thinkstats2.com。

NormalProbabilityPlot is provided by thinkstats2, which is included in the repository for this book. If you are not familiar with normal probability plots, you can read about them in Chapter 5 of Think Stats at http://thinkstats2.com.

图 4-8显示了结果。灰色线条表示对数据进行线性拟合的结果；深色线条表示数据本身。

Figure 4-8 shows the results. The gray lines show a linear model fit to the data; the dark lines show the data.

正态概率图上的直线表明数据来自高斯分布。除了极值处的随机波动外，这些直线基本呈直线，这表明 UG 噪声的频谱就是 UG 噪声。

A straight line on a normal probability plot indicates that the data come from a Gaussian distribution. Except for some random variation at the extremes, these lines are straight, which indicates that the spectrum of UG noise is UG noise.

UU噪声的频谱也近似于UG噪声。事实上，根据中心极限定理，只要分布具有有限的均值和标准差，且样本数量足够大，几乎任何不相关噪声的频谱都近似服从高斯分布。

The spectrum of UU noise is also UG noise, at least approximately. In fact, by the Central Limit Theorem, the spectrum of almost any uncorrelated noise is approximately Gaussian, as long as the distribution has finite mean and standard deviation, and the number of samples is large.

图 4-8.高斯噪声频谱的实部和虚部的正态概率图。

练习

Exercises

这些练习的答案在chap04soln.ipynb……

Solutions to these exercises are in chap04soln.ipynb.

相关性

Correlation

一般来说，变量之间的相关性意味着，如果你知道其中一个变量的值，你就能获得关于另一个变量的一些信息。量化相关性的方法有很多种，但最常用的是皮尔逊积矩相关系数，通常用ρ表示。对于两个变量x和y，每个变量都包含N 个值：

In general, correlation between variables means that if you know the value of one, you have some information about the other. There are several ways to quantify correlation, but the most common is the Pearson product-moment correlation coefficient, usually denoted ρ. For two variables, x and y, that each contain N values:

其中μx和μy_分别是x和y的均值，σx_和σy_{分别是它们}的标准差。

where μ_x and μ_y are the means of x and y, and σ_x and σ_y are their standard deviations.

皮尔逊相关系数 ρ 的值始终介于 -1 和 +1 之间（包含 -1 和 +1）。如果ρ为正，则称相关性为正，这意味着当一个变量较高时，另一个变量也往往较高。如果ρ为负，则相关性为负，这意味着当一个变量较高时，另一个变量往往较低。

Pearson’s correlation is always between –1 and +1 (including both). If ρ is positive, we say that the correlation is positive, which means that when one variable is high, the other tends to be high. If ρ is negative, the correlation is negative, so when one variable is high, the other tends to be low.

ρ的大小表示相关性的强度。如果ρ为 1 或 -1，则变量完全相关，这意味着如果你知道其中一个变量的值，就可以完美地预测另一个变量的值。如果ρ接近于零，则相关性可能很弱，因此如果你知道其中一个变量的值，就无法很好地预测其他变量的值。

The magnitude of ρ indicates the strength of the correlation. If ρ is 1 or –1, the variables are perfectly correlated, which means that if you know one, you can make a perfect prediction about the other. If ρ is near zero, the correlation is probably weak, so if you know one, it doesn’t tell you much about the others.

我说“可能较弱”，是因为也可能存在相关系数无法捕捉到的非线性关系。非线性关系在统计学中通常很重要，但在信号处理中则较少涉及，因此我在此不再赘述。

I say “probably weak” because it is also possible that there is a nonlinear relationship that is not captured by the coefficient of correlation. Nonlinear relationships are often important in statistics, but less often relevant for signal processing, so I won’t say more about them here.

Python 提供了多种计算相关性的方法。np.corrcoef它可以接受任意数量的变量，并计算出一个相关矩阵，其中包含每对变量之间的相关性。

Python provides several ways to compute correlations. np.corrcoef takes any number of variables and computes a correlation matrix that includes correlations between each pair of variables.

我将举一个只有两个变量的例子。首先，我定义一个函数，用于构造具有不同相位偏移的正弦波：

I’ll present an example with only two variables. First, I define a function that constructs sine waves with different phase offsets:

def make_sine(offset):
    signal = thinkdsp.SinSignal(freq=440, offset=offset)
    wave = signal.make_wave(duration=0.5, framerate=10000)
    返回波

def make_sine(offset):
    signal = thinkdsp.SinSignal(freq=440, offset=offset)
    wave = signal.make_wave(duration=0.5, framerate=10000)
    return wave

接下来，我实例化两个具有不同偏移量的波形：

Next I instantiate two waves with different offsets:

wave1 = make_sine(offset=0)
wave2 = make_sine(offset=1)

wave1 = make_sine(offset=0)
wave2 = make_sine(offset=1)

图 5-1显示了这些波的前几个周期的样子。当一个波很高时，另一个波通常也很高，因此我们预期它们是相关的：

Figure 5-1 shows what the first few periods of these waves look like. When one wave is high, the other is usually high, so we expect them to be correlated:

>>> corr_matrix = np.corrcoef(wave1.ys, wave2.ys, ddof=0)
[[ 1. 0.54]
 [ 0.54 1. ]]

>>> corr_matrix = np.corrcoef(wave1.ys, wave2.ys, ddof=0)
[[ 1.    0.54]
 [ 0.54  1.  ]]

该选项ddof=0表示corrcoef应该像上面的等式一样除以N，而不是使用默认值。

The option ddof=0 indicates that corrcoef should divide by N, as in the equation above, rather than use the default, .

结果是一个相关矩阵。第一个元素是wave1与自身的相关系数，始终为 1。类似地，最后一个元素是wave2与自身的相关系数。

The result is a correlation matrix. The first element is the correlation of wave1 with itself, which is always 1. Similarly, the last element is the correlation of wave2 with itself.

非对角线元素包含我们感兴趣的值，即wave1和的相关性wave2。值 0.54 表示相关性强度为中等。

The off-diagonal elements contain the value we’re interested in, the correlation of wave1 and wave2. The value 0.54 indicates that the strength of the correlation is moderate.

图 5-1.两个相位偏移为 1 弧度的正弦波；它们的相相关系数为 0.54。

随着相位偏移增大，相关性逐渐减小，直到波的相位差达到 180 度，此时相关性为 -1。然后相关性增大，直到相位偏移达到 360 度。此时，我们完成了一个完整的循环，相关性为 1。

As the phase offset increases, this correlation decreases until the waves are 180 degrees out of phase, which yields correlation –1. Then it increases until the offset differs by 360 degrees. At that point we have come full circle and the correlation is 1.

图 5-2显示了正弦波的相关性和相位偏移之间的关系。这条曲线的形状应该很熟悉；它是一条余弦曲线。

Figure 5-2 shows the relationship between correlation and phase offset for a sine wave. The shape of that curve should look familiar; it is a cosine.

图 5-2.两个正弦波的相关性与它们之间的相位差的关系。结果是一个余弦函数。

thinkdsp提供一个简单的界面来计算波之间的相关性：

thinkdsp provides a simple interface for computing the correlation between waves:

>>> wave1.corr(wave2) 
0.54

>>> wave1.corr(wave2)
0.54

序列相关性

Serial Correlation

信号通常代表随时间变化的物理量的测量值。例如，我们处理过的声音信号代表电压（或电流）的测量值，这对应于我们感知到的声音的气压变化。

Signals often represent measurements of quantities that vary in time. For example, the sound signals we’ve worked with represent measurements of voltage (or current), which correspond to the changes in air pressure we perceive as sound.

这类测量几乎总是存在序列相关性，即每个元素与其下一个（或前一个）元素之间的相关性。为了计算序列相关性，我们可以对信号进行移位，然后计算移位后的信号与原始信号的相关性：

Measurements like these almost always have serial correlation, which is the correlation between each element and the next (or the previous). To compute serial correlation, we can shift a signal and then compute the correlation of the shifted version with the original:

def serial_corr(wave, lag=1):
    n = len(wave)
    y1 = wave.ys[lag:]
    y2 = wave.ys[:n-lag]
    corr = np.corrcoef(y1, y2, ddof=0)[0, 1]
    返回正确

def serial_corr(wave, lag=1):
    n = len(wave)
    y1 = wave.ys[lag:]
    y2 = wave.ys[:n-lag]
    corr = np.corrcoef(y1, y2, ddof=0)[0, 1]
    return corr

serial_corr函数接受一个Wave对象和lag一个参数，其中参数是要将波形平移的整数位数。它计算波形与其自身平移版本之间的相关性。

serial_corr takes a Wave object and lag, which is the integer number of places to shift the wave. It computes the correlation of the wave with a shifted version of itself.

我们可以用上一章的噪声信号来测试这个函数。根据噪声的生成方式（更不用说它的名称），我们预期 UU 噪声是不相关的：

We can test this function with the noise signals from the previous chapter. We expect UU noise to be uncorrelated, based on the way it’s generated (not to mention the name):

signal = thinkdsp.UncorrelatedGaussianNoise()
wave = signal.make_wave(duration=0.5, framerate=11025)
serial_corr(wave)

signal = thinkdsp.UncorrelatedGaussianNoise()
wave = signal.make_wave(duration=0.5, framerate=11025)
serial_corr(wave)

我运行这个例子后得到 0.006，这表明序列相关性非常小。你运行后可能会得到不同的值，但应该也同样很小。

When I ran this example, I got 0.006, which indicates a very small serial correlation. You might get a different value when you run it, but it should be comparably small.

在布朗噪声信号中，每个值都是前一个值加上一个随机“步长”的总和，因此我们预期会存在很强的序列相关性：

In a Brownian noise signal, each value is the sum of the previous value and a random “step”, so we expect a strong serial correlation:

signal = thinkdsp.BrownianNoise()
wave = signal.make_wave(duration=0.5, framerate=11025)
serial_corr(wave)

signal = thinkdsp.BrownianNoise()
wave = signal.make_wave(duration=0.5, framerate=11025)
serial_corr(wave)

果然，我得到的结果大于 0.999。

Sure enough, the result I got is greater than 0.999.

由于粉红噪声在某种意义上介于布朗噪声和UU噪声之间，我们可以预期它们之间存在中间相关性：

Since pink noise is in some sense between Brownian noise and UU noise, we might expect an intermediate correlation:

signal = thinkdsp.PinkNoise(beta=1)
wave = signal.make_wave(duration=0.5, framerate=11025)
serial_corr(wave)

signal = thinkdsp.PinkNoise(beta=1)
wave = signal.make_wave(duration=0.5, framerate=11025)
serial_corr(wave)

当参数为 时，我得到了序列相关性为 0.851。随着参数从（不相关噪声）变化到（布朗运动），序列相关性从 0 变化到接近 1，如图 5-3所示。

With parameter , I got a serial correlation of 0.851. As we vary the parameter from , which is uncorrelated noise, to , which is Brownian, serial correlation ranges from 0 to almost 1, as shown in Figure 5-3.

图 5-3.粉红噪声的序列相关性，参数范围为 100% 到 100% 不等。

自相关性

Autocorrelation

在上一节中，我们计算了每个值与下一个值之间的相关性，因此我们将数组的元素向右移动了 1。但是，我们可以很容易地计算具有不同滞后的序列相关性。

In the previous section we computed the correlation between each value and the next, so we shifted the elements of the array by 1. But we can easily compute serial correlations with different lags.

你可以将 视为serial_corr一个函数，它将 的每个值映射lag到相应的相关性，我们可以通过遍历 的值来评估该函数lag：

You can think of serial_corr as a function that maps from each value of lag to the corresponding correlation, and we can evaluate that function by looping through values of lag:

def autocorr(wave):
    lags = range(len(wave.ys)//2)
    corrs = [serial_corr(wave, lag) for lag in lags]
    返回滞后，相关系数

def autocorr(wave):
    lags = range(len(wave.ys)//2)
    corrs = [serial_corr(wave, lag) for lag in lags]
    return lags, corrs

autocorr接受一个Wave对象，并返回自相关函数作为一对序列：lags是一个从 0 到波长度一半的整数序列；corrs是每个滞后的序列相关性。

autocorr takes a Wave object and returns the autocorrelation function as a pair of sequences: lags is a sequence of integers from 0 to half the length of the wave; corrs is the sequence of serial correlations for each lag.

图 5-4显示了三种不同β值下粉红噪声的自相关函数。当β值较小时，信号的相关性较弱，自相关函数迅速下降至零。当 β 值较大时，序列相关性更强，下降速度也更慢。即使存在较长的滞后，序列相关性仍然很强；这种现象被称为长程依赖性，因为它表明信号中的每个值都依赖于许多先前的值。

Figure 5-4 shows autocorrelation functions for pink noise with three values of β. For low values of β, the signal is less correlated, and the autocorrelation function drops off to zero quickly. For larger values, serial correlation is stronger and drops off more slowly. With serial correlation is strong even for long lags; this phenomenon is called long-range dependence, because it indicates that each value in the signal depends on many preceding values.

图 5-4.具有一系列参数的粉红噪声的自相关函数。

周期信号的自相关性

Autocorrelation of Periodic Signals

粉红噪声的自相关性具有有趣的数学性质，但应用范围有限。周期信号的自相关性则更有用。

The autocorrelation of pink noise has interesting mathematical properties, but limited applications. The autocorrelation of periodic signals is more useful.

例如，我从https://freesound.org下载了一段鸟鸣声的录音；本书的仓库中包含该文件（https://github.com/AllenDowney/ThinkDSP/blob/master/code/28042__bcjordan__voicedownbew.wav）。您可以使用本章的 Jupyter notebook来播放它。chap05.ipynb

As an example, I downloaded from https://freesound.org a recording of someone singing a chirp; the repository for this book includes the file (https://github.com/AllenDowney/ThinkDSP/blob/master/code/28042__bcjordan__voicedownbew.wav). You can use the Jupyter notebook for this chapter, chap05.ipynb, to play it.

图 5-5显示了该波的频谱图。基频和一些谐波清晰可见。啁啾声起始于 500 Hz 附近，并下降到约 300 Hz，大致位于 C5 到 E4 之间。

Figure 5-5 shows the spectrogram of this wave. The fundamental frequency and some of the harmonics show up clearly. The chirp starts near 500 Hz and drops down to about 300 Hz, roughly from C5 to E4.

为了估算特定时刻的音高，我们可以使用频谱，但效果并不理想。为了说明原因，我将截取波形的一小段并绘制其频谱：

To estimate pitch at a particular point in time, we could use the spectrum, but it doesn’t work very well. To see why not, I’ll take a short segment from the wave and plot its spectrum:

持续时间 = 0.01
segment = wave.segment(start=0.2, duration=duration)
spectrum = segment.make_spectrum()
spectrum.plot(high=1000)

duration = 0.01
segment = wave.segment(start=0.2, duration=duration)
spectrum = segment.make_spectrum()
spectrum.plot(high=1000)

图 5-5.鸣叫声的频谱图。

这段音频片段从 0.2 秒开始，持续 0.01 秒。图 5-6显示了它的频谱。在 400 Hz 附近有一个明显的峰值，但很难精确地确定音高。该片段的长度为 441 个采样点，帧率为 44,100 Hz，因此频率分辨率为 100 Hz（参见“Gabor 极限”）。这意味着估计的音高可能存在 50 Hz 的偏差；在音乐术语中，350 Hz 到 450 Hz 的范围大约是 5 个半音，这是一个很大的差异！

This segment starts at 0.2 seconds and lasts 0.01 seconds. Figure 5-6 shows its spectrum. There is a clear peak near 400 Hz, but it is hard to identify the pitch precisely. The length of the segment is 441 samples at a frame rate of 44,100 Hz, so the frequency resolution is 100 Hz (see “The Gabor Limit”). That means the estimated pitch might be off by 50 Hz; in musical terms, the range from 350 Hz to 450 Hz is about 5 semitones, which is a big difference!

图 5-6.人声啁啾声片段的频谱。

我们可以通过取更长的片段来获得更好的频率分辨率，但由于音调随时间变化，我们也会出现“运动模糊”；也就是说，峰值会在片段的开始音调和结束音调之间扩散，正如我们在“啁啾声的频谱”中所看到的那样。

We could get better frequency resolution by taking a longer segment, but since the pitch is changing over time, we would also get “motion blur”; that is, the peak would spread between the start and end pitch of the segment, as we saw in “Spectrum of a Chirp”.

我们可以使用自相关来更精确地估计音高。如果信号是周期性的，我们预期当滞后时间等于周期时，自相关函数会出现峰值。

We can estimate pitch more precisely using autocorrelation. If a signal is periodic, we expect the autocorrelation to spike when the lag equals the period.

为了说明这种方法的有效性，我将绘制同一段录音中的两个片段：

To show why that works, I’ll plot two segments from the same recording:

def plot_shifted(wave, offset=0.001, start=0.2):
    thinkplot.preplot(2)
    segment1 = wave.segment(start=start, duration=0.01)
    segment1.plot(linewidth=2, alpha=0.8)

    segment2 = wave.segment(start=start-offset, duration=0.01)
    segment2.shift(offset)
    segment2.plot(linewidth=2, alpha=0.4)

    corr = segment1.corr(segment2)
    text = r'$\rho =$ %.2g' % corr
    thinkplot.text(segment1.start+0.0005, -0.8, text)
    thinkplot.config(xlabel='时间（秒）')

def plot_shifted(wave, offset=0.001, start=0.2):
    thinkplot.preplot(2)
    segment1 = wave.segment(start=start, duration=0.01)
    segment1.plot(linewidth=2, alpha=0.8)

    segment2 = wave.segment(start=start-offset, duration=0.01)
    segment2.shift(offset)
    segment2.plot(linewidth=2, alpha=0.4)

    corr = segment1.corr(segment2)
    text = r'$\rho =$ %.2g' % corr
    thinkplot.text(segment1.start+0.0005, -0.8, text)
    thinkplot.config(xlabel='Time (s)')

一段信号起始于 0.2 秒；另一段信号起始于前者 0.0023 秒之后。图 5-7显示了结果。这两段信号非常相似，相关性为 0.99。该结果表明周期接近 0.0023 秒，对应的频率为 435 Hz。

One segment starts at 0.2 seconds; the other starts 0.0023 seconds later. Figure 5-7 shows the result. The segments are similar, and their correlation is 0.99. This result suggests that the period is near 0.0023 seconds, which corresponds to a frequency of 435 Hz.

图 5-7.啁啾声的两个片段，一个片段比另一个片段晚 0.0023 秒开始。

在这个例子中，我通过反复试验来估算周期。为了实现自动化，我们可以使用自相关函数：

For this example, I estimated the period by trial and error. To automate the process, we can use the autocorrelation function:

滞后，相关系数 = autocorr(segment)
thinkplot.plot(lags, corrs)

lags, corrs = autocorr(segment)
thinkplot.plot(lags, corrs)

图 5-8显示了从秒开始的片段的自相关函数。第一个峰值出现在lag=101。我们可以按如下方式计算与该周期对应的频率：

Figure 5-8 shows the autocorrelation function for the segment starting at seconds. The first peak occurs at lag=101. We can compute the frequency that corresponds to that period like this:

周期 = 延迟 / 片段帧率
频率 = 1 / 周期

period = lag / segment.framerate
frequency = 1 / period

估计的基频为 437 Hz。为了评估估计的精度，我们可以分别使用 100 和 102 的滞后进行相同的计算，这对应于 432 Hz 和 441 Hz 的频率。使用自相关法得到的频率精度小于 10 Hz，而使用频谱法得到的精度为 100 Hz。用音乐术语来说，预期误差约为 30 音分（三分之一半音）。

The estimated fundamental frequency is 437 Hz. To evaluate the precision of the estimate, we can run the same computation with lags 100 and 102, which correspond to frequencies 432 and 441 Hz. The frequency precision using autocorrelation is less than 10 Hz, compared with 100 Hz using the spectrum. In musical terms, the expected error is about 30 cents (a third of a semitone).

图 5-8.啁啾信号片段的自相关函数。

相关性即点积

Correlation as Dot Product

本章开头，我给出了皮尔逊相关系数的定义：

I started the chapter with this definition of Pearson’s correlation coefficient:

然后我用ρ定义了序列相关和自相关。这与统计学中这些术语的用法一致，但在信号处理的背景下，它们的定义略有不同。

Then I used ρ to define serial correlation and autocorrelation. That’s consistent with how these terms are used in statistics, but in the context of signal processing, the definitions are a little different.

在信号处理中，我们经常处理均值为 0 的无偏信号和标准差为 1 的归一化信号。在这种情况下，ρ的定义简化为：

In signal processing, we are often working with unbiased signals, where the mean is 0, and normalized signals, where the standard deviation is 1. In that case, the definition of ρ simplifies to:

而且通常还会进一步简化：

And it is common to simplify even further:

这种相关性的定义并非“标准化”的，因此其值通常不在 -1 到 1 之间。但它具有其他有用的特性。

This definition of correlation is not “standardized”, so it doesn’t generally fall between –1 and 1. But it has other useful properties.

如果你把x和y看作向量，你可能会认出这个公式是点积公式。参见http://en.wikipedia.org/wiki/Dot_product。

If you think of x and y as vectors, you might recognize this formula as the dot product, . See http://en.wikipedia.org/wiki/Dot_product.

点积表示信号的相似程度。如果将它们归一化，使其标准差为 1：

The dot product indicates the degree to which the signals are similar. If they are normalized so their standard deviations are 1:

其中θ是向量之间的夹角。这就解释了为什么图5-2是一条余弦曲线。

where θ is the angle between the vectors. And that explains why Figure 5-2 is a cosine curve.

使用 NumPy

Using NumPy

NumPy 提供了一个函数 ` correlateccorrelation`，用于计算两个函数的相关性或一个函数的自相关性。我们可以用它来计算上一节中片段的自相关性：

NumPy provides a function, correlate, that computes the correlation of two functions or the autocorrelation of one function. We can use it to compute the autocorrelation of the segment from the previous section:

corrs2 = np.correlate(segment.ys, segment.ys, mode='same')

corrs2 = np.correlate(segment.ys, segment.ys, mode='same')

该选项mode指定correlate要使用的范围lag。当值为 时'same'，范围为到，其中N为波阵列的长度。

The option mode tells correlate what range of lag to use. With the value 'same', the range is from to , where N is the length of the wave array.

图 5-9显示了结果。由于两个信号完全相同，因此它是对称的，一个信号的负滞后与另一个信号的正滞后具有相同的效果。为了与之前的结果进行比较autocorr，我们可以选择后半部分：

Figure 5-9 shows the result. It is symmetric because the two signals are identical, so a negative lag on one has the same effect as a positive lag on the other. To compare with the results from autocorr, we can select the second half:

N = len(corrs2)
half = corrs2[N//2:]

N = len(corrs2)
half = corrs2[N//2:]

图 5-9.使用 np.correlate 计算的自相关函数。

比较图 5-9和图 5-8，你会发现，随着滞后阶数的增加，计算得到的关联度np.correlate逐渐减小。这是因为该方法np.correlate使用了非标准化的关联度定义；随着滞后阶数的增大，两个信号重叠区域的点数减少，因此关联度的幅度也随之减小。

If you compare Figure 5-9 to Figure 5-8, you’ll notice that the correlations computed by np.correlate get smaller as the lags increase. That’s because np.correlate uses the unstandardized definition of correlation; as the lag gets bigger, the number of points in the overlap between the two signals gets smaller, so the magnitude of the correlations decreases.

我们可以通过除以长度来修正这个问题：

We can correct that by dividing through by the lengths:

长度 = range(N, N//2, -1)
一半/=长度

lengths = range(N, N//2, -1)
half /= lengths

最后，我们可以对结果进行归一化处理，使相关性为lag=01：

Finally, we can normalize the results so the correlation with lag=0 is 1:

half /= half[0]

half /= half[0]

autocorr经过这些调整，两种方法计算的结果np.correlate几乎相同，仅有1%到2%的差异。原因并不重要，但如果您好奇的话：autocorr对每个滞后阶数的相关性进行独立标准化；而对于np.correlate，我们最后对所有滞后阶数进行了标准化。

With these adjustments, the results computed by autocorr and np.correlate are nearly the same. They still differ by 1–2%. The reason is not important, but if you are curious: autocorr standardizes the correlations independently for each lag; for np.correlate, we standardized them all at the end.

更重要的是，现在你知道什么是自相关，如何使用它来估计信号的基本周期，以及计算自相关的两种方法。

More importantly, now you know what autocorrelation is, how to use it to estimate the fundamental period of a signal, and two ways to compute it.

合成

Synthesis

假设我给你一个振幅列表和一个频率列表，并要求你构建一个信号，该信号是这些频率分量的总和。使用thinkdsp模块中的对象，有一种简单的方法可以执行此操作，这称为合成：

Suppose I give you a list of amplitudes and a list of frequencies, and ask you to construct a signal that is the sum of these frequency components. Using objects in the thinkdsp module, there is a simple way to perform this operation, which is called synthesis:

def synthesize1(amps, fs, ts):
    components = [thinkdsp.CosSignal(freq, amp)
                  对于 zip(amps, fs) 中的每个 amp 和 freq]
    signal = thinkdsp.SumSignal(*components)

    ys = signal.evaluate(ts)
    返回 ys

def synthesize1(amps, fs, ts):
    components = [thinkdsp.CosSignal(freq, amp)
                  for amp, freq in zip(amps, fs)]
    signal = thinkdsp.SumSignal(*components)

    ys = signal.evaluate(ts)
    return ys

amps是振幅列表，fs是频率列表，ts是信号评估时间序列。

amps is a list of amplitudes, fs is the list of frequencies, and ts is the sequence of times where the signal should be evaluated.

components是一个对象列表CosSignal，每个对象对应一个振幅-频率对。SumSignal表示这些频率分量的总和。

components is a list of CosSignal objects, one for each amplitude–frequency pair. SumSignal represents the sum of these frequency components.

最后，evaluate计算每个时刻的信号值ts。

Finally, evaluate computes the value of the signal at each time in ts.

我们可以这样测试这个函数：

We can test this function like this:

amps = np.array([0.6, 0.25, 0.1, 0.05])
fs = [100, 200, 300, 400]
帧率 = 11025

ts = np.linspace(0, 1, framerate)
ys = synthesize1(amps, fs, ts)
wave = thinkdsp.Wave(ys, framerate)

amps = np.array([0.6, 0.25, 0.1, 0.05])
fs = [100, 200, 300, 400]
framerate = 11025

ts = np.linspace(0, 1, framerate)
ys = synthesize1(amps, fs, ts)
wave = thinkdsp.Wave(ys, framerate)

此示例生成一个包含 100 Hz 基频和三个谐波的信号（100 Hz 是一个尖锐的 G2 谐波）。它以每秒 11,025 帧的速度渲染该信号 1 秒钟，并将结果放入一个Wave对象中。

This example makes a signal that contains a fundamental frequency at 100 Hz and three harmonics (100 Hz is a sharp G2). It renders the signal for 1 second at 11,025 frames per second and puts the results into a Wave object.

从概念上讲，合成非常简单。但这种形式对分析帮助不大，分析是它的逆问题：给定一个波，我们如何识别其频率分量及其振幅？

Conceptually, synthesis is pretty simple. But in this form it doesn’t help much with analysis, which is the inverse problem: given the wave, how do we identify the frequency components and their amplitudes?

基于阵列的合成

Synthesis with Arrays

另一种写法是synthesize：

Here’s another way to write synthesize:

def synthesize2(amps, fs, ts):
    args = np.outer(ts, fs)
    M = np.cos(PI2 * args)
    ys = np.dot(M, amps)
    返回 ys

def synthesize2(amps, fs, ts):
    args = np.outer(ts, fs)
    M = np.cos(PI2 * args)
    ys = np.dot(M, amps)
    return ys

这个函数看起来完全不同，但功能相同。我们来看看它是如何工作的：

This function looks very different, but it does the same thing. Let’s see how it works:

1. np.outer计算和的外积。结果是一个数组，其中每一行对应 的一个元素，每一列对应 的一个元素。数组中的每个元素都是频率和时间的乘积。tsfstsfs

2. np.outer computes the outer product of ts and fs. The result is an array with one row for each element of ts and one column for each element of fs. Each element in the array is the product of a frequency and a time, .

3. 我们乘以args并应用cos，因此结果的每个元素都是。由于沿着ts列向下运行，因此每一列都包含一个特定频率的余弦信号，该信号在一系列时间点上进行计算。

4. We multiply args by and apply cos, so each element of the result is . Since the ts runs down the columns, each column contains a cosine signal at a particular frequency, evaluated at a sequence of times.

5. np.dot将矩阵的每一行逐元素地M乘以amps向量，然后将乘积相加。用线性代数的术语来说，我们是在将矩阵M乘以向量。用信号学的术语来说，我们是在计算频率分量amps的加权和。

6. np.dot multiplies each row of M by amps, elementwise, and then adds up the products. In terms of linear algebra, we are multiplying a matrix, M, by a vector, amps. In terms of signals, we are computing the weighted sum of frequency components.

图 6-1显示了该计算的结构。矩阵的每一行M对应于 0.0 到 1.0 秒之间的时间；t _{_n}是第n行的时间。每一列对应于 100 到 400 Hz 之间的频率；f _{_k}是第k列的频率。

Figure 6-1 shows the structure of this computation. Each row of the matrix, M, corresponds to a time from 0.0 to 1.0 seconds; t_n is the time of the nth row. Each column corresponds to a frequency from 100 to 400 Hz; f_k is the frequency of the kth column.

图 6-1.阵列合成。

我用字母a到d标记了第n行；例如，a的值为。

I labeled the nth row with the letters a through d; as an example, the value of a is .

点积的结果ys是一个向量，其中每个元素对应于矩阵的每一行M。第n个元素，标记为e，是各乘积之和：

The result of the dot product, ys, is a vector with one element for each row of M. The nth element, labeled e, is the sum of products:

同样的道理也适用于其他元素ys。因此，的每个元素ys都是四个频率分量在某一时刻的取值之和，并乘以相应的振幅。这正是我们想要的。

And likewise with the other elements of ys. So each element of ys is the sum of four frequency components, evaluated at a point in time and multiplied by the corresponding amplitudes. And that’s exactly what we wanted.

我们可以使用上一节中的代码来验证这两个版本是否synthesize产生相同的结果：

We can use the code from the previous section to check that the two versions of synthesize produce the same results:

ys1 = synthesize1(amps, fs, ts)
ys2 = synthesize2(amps, fs, ts)
max(abs(ys1 - ys2))

ys1 = synthesize1(amps, fs, ts)
ys2 = synthesize2(amps, fs, ts)
max(abs(ys1 - ys2))

ys1和之间最大的区别ys2约为1e-13，这是我们由于浮点误差而预期的结果。

The biggest difference between ys1 and ys2 is about 1e-13, which is what we expect due to floating-point errors.

用线性代数来编写这个计算可以简化代码，提高运行速度。线性代数为矩阵和向量的运算提供了简洁的表示法。例如，我们可以synthesize这样写：

Writing this computation in terms of linear algebra makes the code smaller and faster. Linear algebra provides concise notation for operations on matrices and vectors. For example, we could write synthesize like this:

其中a是振幅向量，t是时间向量，f是频率向量，⊗是两个向量外积的符号。

where a is a vector of amplitudes, t is a vector of times, f is a vector of frequencies, and ⊗ is the symbol for the outer product of two vectors.

分析

Analysis

现在我们准备解决分析问题。假设我给你一个波，并告诉你它是给定频率集合的余弦波之和。你如何求出每个频率分量的振幅？换句话说，给定频率ys、ts频率 和fs频率 ，你能恢复出频率amps吗？

Now we are ready to solve the analysis problem. Suppose I give you a wave and tell you that it is the sum of cosines with a given set of frequencies. How would you find the amplitude for each frequency component? In other words, given ys, ts, and fs, can you recover amps?

就线性代数而言，第一步与综合相同：我们计算。然后我们需要找到一个，使得；换句话说，我们需要求解一个线性方程组。NumPy 提供了函数linalg.solve，它正是用来做这件事的。

In terms of linear algebra, the first step is the same as for synthesis: we compute . Then we want to find a so that ; in other words, we want to solve a linear system. NumPy provides linalg.solve, which does exactly that.

以下是代码：

Here’s what the code looks like:

def analyze1(ys, fs, ts):
    args = np.outer(ts, fs)
    M = np.cos(PI2 * args)
    amps = np.linalg.solve(M, ys)
    回流放大器

def analyze1(ys, fs, ts):
    args = np.outer(ts, fs)
    M = np.cos(PI2 * args)
    amps = np.linalg.solve(M, ys)
    return amps

前两行使用ts和fs来构建矩阵，M然后np.linalg.solve计算amps。

The first two lines use ts and fs to build the matrix, M. Then np.linalg.solve computes amps.

但这里有个问题。一般来说，只有当矩阵是方阵时，我们才能解线性方程组；也就是说，只有当方程的数量（行数）与未知数的数量（列数）相同时，我们才能解线性方程组。

But there’s a hitch. In general we can only solve a system of linear equations if the matrix is square; that is, if the number of equations (rows) is the same as the number of unknowns (columns).

在这个例子中，我们只有 4 个频率，但我们对信号进行了 11025 次评估。因此，我们的方程数量远远多于未知数的数量。

In this example, we have only 4 frequencies, but we evaluated the signal at 11,025 times. So we have many more equations than unknowns.

一般来说，如果ys包含超过 4 个元素，我们不太可能仅使用 4 个频率对其进行分析。

In general if ys contains more than 4 elements, it is unlikely that we can analyze it using only 4 frequencies.

但在这种情况下，我们知道这些ys值实际上是通过仅添加 4 个频率分量生成的，因此我们可以使用波阵列中的任意 4 个值来恢复amps。

But in this case, we know that the ys values were actually generated by adding only 4 frequency components, so we can use any 4 values from the wave array to recover amps.

为简单起见，我将使用信号的前 4 个样本。利用上一节中得到的 、 和 的值，我们ys可以fs这样运行：tsanalyze1

For simplicity, I’ll use the first 4 samples from the signal. Using the values of ys, fs, and ts from the previous section, we can run analyze1 like this:

n = len(fs)
amps2 = analyze1(ys[:n], fs, ts[:n])

n = len(fs)
amps2 = analyze1(ys[:n], fs, ts[:n])

果不其然，amps2正是如此：

And sure enough, amps2 is:

[ 0.6 0.25 0.1 0.05 ]

[ 0.6   0.25  0.1   0.05 ]

这个算法虽然可行，但速度很慢。求解线性方程组所需的时间与n³^成正比，其中n是矩阵M的列数。我们可以做得更好。

This algorithm works, but it is slow. Solving a linear system of equations takes time proportional to n³, where n is the number of columns in M. We can do better.

正交矩阵

Orthogonal Matrices

求解线性方程组的一种方法是求矩阵的逆。方阵M的逆记为 I ，其性质为 I = 1。I是单位矩阵，其对角线元素均为 1，其余元素均为 0。

One way to solve linear systems is by inverting matrices. The inverse of a square matrix M is written , and it has the property that . I is the identity matrix, which has the value 1 on all diagonal elements and 0 everywhere else.

因此，为了解这个方程，我们可以将等式两边同时乘以，得到：

So, to solve the equation , we can multiply both sides by , which yields:

右侧我们可以用I替换：

On the right side, we can replace with I:

如果我们将I乘以任意向量a，结果仍然是a，所以：

If we multiply I by any vector a, the result is a, so:

这意味着，如果我们能够高效地进行计算，我们就可以通过简单的矩阵乘法（使用）找到一个np.dot。这需要与n²^成正比的时间，比n³更优。

This implies that if we can compute efficiently, we can find a with a simple matrix multiplication (using np.dot). That takes time proportional to n², which is better than n³.

通常来说，矩阵求逆速度较慢，但​​某些特殊情况下速度会更快。特别是，如果矩阵M是正交矩阵，则M的逆矩阵就是M的转置矩阵，记作M ^T。在 NumPy 中，数组转置是一个常数时间操作。它实际上并不移动数组元素；而是创建一个“视图”，改变元素的访问方式。

Inverting a matrix is slow, in general, but some special cases are faster. In particular, if M is orthogonal, the inverse of M is just the transpose of M, written M^T. In NumPy transposing an array is a constant-time operation. It doesn’t actually move the elements of the array; instead, it creates a “view” that changes the way the elements are accessed.

再次强调，如果一个矩阵的转置也是它的逆矩阵，则称该矩阵为正交矩阵；也就是说，。这意味着，我们可以通过计算M ^T M来检查一个矩阵是否正交。

Again, a matrix is orthogonal if its transpose is also its inverse; that is, . That implies that , which means we can check whether a matrix is orthogonal by computing M^TM.

那么让我们看看矩阵在 中的样子synthesize2。在前面的例子中，矩阵 M有 11,025 行，所以最好使用一个更小的例子：

So let’s see what the matrix looks like in synthesize2. In the previous example, M has 11,025 rows, so it might be a good idea to work with a smaller example:

def test1():
    amps = np.array([0.6, 0.25, 0.1, 0.05])
    N = 4.0
    时间单位 = 0.001
    ts = np.arange(N) / N * time_unit
    最大频率 = N / 时间单位 / 2
    fs = np.arange(N) / N * max_freq
    ys = synthesize2(amps, fs, ts)

def test1():
    amps = np.array([0.6, 0.25, 0.1, 0.05])
    N = 4.0
    time_unit = 0.001
    ts = np.arange(N) / N * time_unit
    max_freq = N / time_unit / 2
    fs = np.arange(N) / N * max_freq
    ys = synthesize2(amps, fs, ts)

amps是之前看到的同一幅值向量。由于我们有 4 个频率分量，我们将对信号进行 4 次采样。这样，M就是平方。

amps is the same vector of amplitudes we saw before. Since we have 4 frequency components, we’ll sample the signal at 4 points in time. That way, M is square.

ts是一个向量，其中包含从 0 到 1 个时间单位的等间隔采样时间。我选择的时间单位是 1 毫秒，但这只是任意选择，稍后我们会看到，它最终会被排除在计算之外。

ts is a vector of equally spaced sample times in the range from 0 to 1 time unit. I chose the time unit to be 1 millisecond, but it is an arbitrary choice, and we will see in a minute that it drops out of the computation anyway.

由于帧速率为每时间单位N 个样本，因此奈奎斯特频率为N / time_unit / 2，在本例中为 2000 Hz。所以fs是一个介于 0 和 2000 Hz 之间的等间隔频率向量。

Since the frame rate is N samples per time unit, the Nyquist frequency is N / time_unit / 2, which is 2000 Hz in this example. So fs is a vector of equally spaced frequencies between 0 and 2000 Hz.

使用这些值ts，fs矩阵M为：

With these values of ts and fs, the matrix, M, is:

[[ 1. 1. 1. 1. ]
 [ 1. 0.707 0. -0.707]
 [ 1. 0. -1. -0. ]
 [ 1. -0.707 -0. 0.707]

[[ 1.     1.     1.     1.   ]
 [ 1.     0.707  0.    -0.707]
 [ 1.     0.    -1.    -0.   ]
 [ 1.    -0.707 -0.     0.707]]

你可能认出 0.707 是 的近似值，即。你可能还会注意到这个矩阵是对称的，这意味着矩阵中第 i 个元素始终等于第 j 个元素。这意味着M是它自身的转置；也就是说，。

You might recognize 0.707 as an approximation of , which is . You also might notice that this matrix is symmetric, which means that the element at always equals the element at . This implies that M is its own transpose; that is, .

但遗憾的是，M不是正交的。如果我们计算M ^T M，我们会得到：

But sadly, M is not orthogonal. If we compute M^TM, we get:

[[ 4. 1. -0. 1.]
 [ 1. 2. 1. -0.]
 [-0. 1. 2. 1.]
 [ 1. -0. 1. 2.]

[[ 4.  1. -0.  1.]
 [ 1.  2.  1. -0.]
 [-0.  1.  2.  1.]
 [ 1. -0.  1.  2.]]

那不是单位矩阵。

And that’s not the identity matrix.

DCT-IV

DCT-IV

但是，如果我们选择ts并fs仔细考虑，我们可以使M正交。有几种方法可以做到这一点，这也是为什么离散余弦变换 (DCT) 有多个版本的原因。

But if we choose ts and fs carefully, we can make M orthogonal. There are several ways to do it, which is why there are several versions of the Discrete Cosine Transform (DCT).

一个简单的办法是将变速器移位ts半个fs单位。这个版本叫做DCT-IV，其中“IV”是罗马数字，表示这是DCT变速器八个版本中的第四个。

One simple option is to shift ts and fs by a half unit. This version is called DCT-IV, where “IV” is a roman numeral indicating that this is the fourth of eight versions of the DCT.

以下是更新版本test1：

Here’s an updated version of test1:

def test2():
    amps = np.array([0.6, 0.25, 0.1, 0.05])
    N = 4.0
    ts = (0.5 + np.arange(N)) / N
    fs = (0.5 + np.arange(N)) / 2
    ys = synthesize2(amps, fs, ts)

def test2():
    amps = np.array([0.6, 0.25, 0.1, 0.05])
    N = 4.0
    ts = (0.5 + np.arange(N)) / N
    fs = (0.5 + np.arange(N)) / 2
    ys = synthesize2(amps, fs, ts)

如果将此版本与之前的版本进行比较，你会注意到两处变化。首先，我给ts和 都加上了 0.5 fs。其次，我约掉了time_units，这简化了 的表达式fs。

If you compare this to the previous version, you’ll notice two changes. First, I added 0.5 to ts and fs. Second, I canceled out time_units, which simplifies the expression for fs.

根据这些值，M为：

With these values, M is:

[[ 0.981 0.831 0.556 0.195]
 [ 0.831 -0.195 -0.981 -0.556]
 [ 0.556 -0.981 0.195 0.831]
 [ 0.195 -0.556 0.831 -0.981]

[[ 0.981  0.831  0.556  0.195]
 [ 0.831 -0.195 -0.981 -0.556]
 [ 0.556 -0.981  0.195  0.831]
 [ 0.195 -0.556  0.831 -0.981]]

而M ^T M是：

And M^TM is:

[[ 2. 0. 0. 0.]
 [ 0. 2. -0. 0.]
 [ 0. -0. 2. -0.]
 [ 0. 0. -0. 2.]

[[ 2.  0.  0.  0.]
 [ 0.  2. -0.  0.]
 [ 0. -0.  2. -0.]
 [ 0.  0. -0.  2.]]

部分非对角线元素显示为 -0，这意味着其浮点表示是一个较小的负数。该矩阵非常接近2I，这意味着M几乎是正交的；它仅相差 2 倍。就我们的目的而言，这已经足够了。

Some of the off-diagonal elements are displayed as –0, which means that the floating-point representation is a small negative number. This matrix is very close to 2I, which means M is almost orthogonal; it’s just off by a factor of 2. And for our purposes, that’s good enough.

因为M是对称的且（几乎）正交的，所以M的逆矩阵就是 M' 。现在我们可以写出一个更高效的版本analyze：

Because M is symmetric and (almost) orthogonal, the inverse of M is just . Now we can write a more efficient version of analyze:

def analyze2(ys, fs, ts):
    args = np.outer(ts, fs)
    M = np.cos(PI2 * args)
    amps = np.dot(M, ys) / 2
    回流放大器

def analyze2(ys, fs, ts):
    args = np.outer(ts, fs)
    M = np.cos(PI2 * args)
    amps = np.dot(M, ys) / 2
    return amps

我们不用np.linalg.solve，而是直接乘以。

Instead of using np.linalg.solve, we just multiply by .

结合test2以上两点analyze2，我们可以编写 DCT-IV 的实现：

Combining test2 and analyze2, we can write an implementation of DCT-IV:

def dct_iv(ys):
    N = len(ys)
    ts = (0.5 + np.arange(N)) / N
    fs = (0.5 + np.arange(N)) / 2
    args = np.outer(ts, fs)
    M = np.cos(PI2 * args)
    amps = np.dot(M, ys) / 2
    回流放大器

def dct_iv(ys):
    N = len(ys)
    ts = (0.5 + np.arange(N)) / N
    fs = (0.5 + np.arange(N)) / 2
    args = np.outer(ts, fs)
    M = np.cos(PI2 * args)
    amps = np.dot(M, ys) / 2
    return amps

再次强调，ys这是波阵列。我们不必将ts和fs作为参数传递；dct_iv可以根据N和 的长度来计算它们ys。

Again, ys is the wave array. We don’t have to pass ts and fs as parameters; dct_iv can figure them out based on N, the length of ys.

如果我们的理解正确，这个函数应该能够解决分析问题；也就是说，给定一个条件，ys它应该能够恢复结果amps。我们可以这样测试它：

If we’ve got it right, this function should solve the analysis problem; that is, given ys it should be able to recover amps. We can test it like this:

amps = np.array([0.6, 0.25, 0.1, 0.05])
N = 4.0
ts = (0.5 + np.arange(N)) / N
fs = (0.5 + np.arange(N)) / 2
ys = synthesize2(amps, fs, ts)
amps2 = dct_iv(ys)
max(abs(amps - amps2))

amps = np.array([0.6, 0.25, 0.1, 0.05])
N = 4.0
ts = (0.5 + np.arange(N)) / N
fs = (0.5 + np.arange(N)) / 2
ys = synthesize2(amps, fs, ts)
amps2 = dct_iv(ys)
max(abs(amps - amps2))

首先amps，我们合成一个波阵列，然后用它dct_iv来计算。和amps2之间的最大差异约为，这符合我们预期，因为存在浮点误差。ampsamps21e-16

Starting with amps, we synthesize a wave array, then use dct_iv to compute amps2. The biggest difference between amps and amps2 is about 1e-16, which is what we expect due to floating-point errors.

逆离散余弦变换

Inverse DCT

最后，注意到analyze2和synthesize2几乎完全相同。唯一的区别在于analyze2将结果除以 2。我们可以利用这一特性来计算逆 DCT：

Finally, notice that analyze2 and synthesize2 are almost identical. The only difference is that analyze2 divides the result by 2. We can use this insight to compute the inverse DCT:

def inverse_dct_iv(amps):
    返回 dct_iv(amps) * 2

def inverse_dct_iv(amps):
    return dct_iv(amps) * 2

inverse_dct_iv它解决了合成问题：它接受振幅向量并返回波数组ys。我们可以通过从 开始amps，应用inverse_dct_iv和dct_iv，并测试是否能得到与初始值相同的结果来验证它：

inverse_dct_iv solves the synthesis problem: it takes the vector of amplitudes and returns the wave array, ys. We can test it by starting with amps, applying inverse_dct_iv and dct_iv, and testing that we get back what we started with:

安培数 = [0.6, 0.25, 0.1, 0.05]
ys = inverse_dct_iv(amps)
amps2 = dct_iv(ys)
max(abs(amps - amps2))

amps = [0.6, 0.25, 0.1, 0.05]
ys = inverse_dct_iv(amps)
amps2 = dct_iv(ys)
max(abs(amps - amps2))

再次强调，最大的区别在于1e-16……

Again, the biggest difference is about 1e-16.

DCT 类

The Dct Class

thinkdsp提供了一个类，它以与封装 FFTDct相同的方式封装了 DCT 。要创建一个对象，您可以调用以下方法：SpectrumDctmake_dctWave

thinkdsp provides a Dct class that encapsulates the DCT in the same way the Spectrum class encapsulates the FFT. To make a Dct object, you can invoke make_dct on a Wave:

signal = thinkdsp.TriangleSignal(freq=400)
wave = signal.make_wave(duration=1.0, framerate=10000)
dct = wave.make_dct()
dct.plot()

signal = thinkdsp.TriangleSignal(freq=400)
wave = signal.make_wave(duration=1.0, framerate=10000)
dct = wave.make_dct()
dct.plot()

结果为 400 Hz 三角波的 DCT，如图 6-2所示。DCT 的值可以是正值或负值；DCT 中的负值对应于取反余弦，或者等价地，对应于移位 180 度的余弦。

The result is the DCT of a triangle wave at 400 Hz, shown in Figure 6-2. The values of the DCT can be positive or negative; a negative value in the DCT corresponds to a negated cosine or, equivalently, to a cosine shifted by 180 degrees.

图 6-2. 400 Hz 三角波信号的 DCT，采样率为 10 kHz。

make_dct采用 DCT-II，这是最常见的 DCT 类型，由以下公司提供scipy.fftpack：

make_dct uses DCT-II, which is the most common type of DCT, provided by scipy.fftpack:

导入 scipy.fftpack

# Wave 类：
    def make_dct(self):
        N = len(self.ys)
        hs = scipy.fftpack.dct(self.ys, type=2)
        fs = (0.5 + np.arange(N)) / 2
        返回 Dct(hs, fs, self.framerate)

import scipy.fftpack

# class Wave:
    def make_dct(self):
        N = len(self.ys)
        hs = scipy.fftpack.dct(self.ys, type=2)
        fs = (0.5 + np.arange(N)) / 2
        return Dct(hs, fs, self.framerate)

计算结果dct存储在 中hs。按照“DCT-IV”方法计算的相应频率存储在 中fs。然后，两者均用于初始化Dct对象。

The results from dct are stored in hs. The corresponding frequencies, computed as in “DCT-IV”, are stored in fs. And then both are used to initialize the Dct object.

Dct提供make_wave执行逆DCT的函数。我们可以这样测试它：

Dct provides make_wave, which performs the inverse DCT. We can test it like this:

wave2 = dct.make_wave()
max(abs(wave.ys-wave2.ys))

wave2 = dct.make_wave()
max(abs(wave.ys-wave2.ys))

ys1和之间最大的区别ys2约为1e-16，这是我们由于浮点误差而预期的结果。

The biggest difference between ys1 and ys2 is about 1e-16, which is what we expect due to floating-point errors.

make_wave用途scipy.fftpack.idct：

make_wave uses scipy.fftpack.idct:

# 类 Dct
    def make_wave(self):
        n = len(self.hs)
        ys = scipy.fftpack.idct(self.hs, type=2) / 2 / n
        返回 Wave(ys, framerate=self.framerate)

# class Dct
    def make_wave(self):
        n = len(self.hs)
        ys = scipy.fftpack.idct(self.hs, type=2) / 2 / n
        return Wave(ys, framerate=self.framerate)

默认情况下，逆 DCT 不会对结果进行归一化，因此我们必须除以2N。

By default, the inverse DCT doesn’t normalize the result, so we have to divide through by 2N.

练习

Exercises

以下练习中，我提供了一些初始代码chap06starter.ipynb。答案在chap06soln.ipynb。

For the following exercises, I provide some starter code in chap06starter.ipynb. Solutions are in chap06soln.ipynb.

复指数

Complex Exponentials

数学中最有趣的操作之一是将一种运算从一种类型推广到另一种类型。例如，阶乘是一个作用于整数的函数；n的阶乘的自然定义是从 1 到n 的所有整数的乘积。

One of the more interesting moves in mathematics is the generalization of an operation from one type to another. For example, a factorial is a function that operates on integers; the natural definition for the factorial of n is the product of all integers from 1 to n.

如果你对计算阶乘感兴趣，你可能会想知道如何计算像 3.5 这样的非整数的阶乘。由于自然定义不适用，你可能会寻找其他方法来计算阶乘函数，即适用于非整数的方法。

If you are of a certain inclination, you might wonder how to compute the factorial of a non-integer like 3.5. Since the natural definition doesn’t apply, you might look for other ways to compute the factorial function, ways that would work with non-integers.

1730 年，莱昂哈德·欧拉发现了一个函数，它是阶乘函数的推广，我们称之为伽玛函数（参见http://en.wikipedia.org/wiki/Gamma_function）。

In 1730, Leonhard Euler found one, a generalization of the factorial function that we know as the gamma function (see http://en.wikipedia.org/wiki/Gamma_function).

欧拉还发现了应用数学中最有用的推广之一——复指数函数。

Euler also found one of the most useful generalizations in applied mathematics, the complex exponential function.

幂运算的自然定义是重复乘法；例如，。但这个定义不适用于非整数指数。

The natural definition of exponentiation is repeated multiplication; for example, . But this definition doesn’t apply to non-integer exponents.

然而，指数运算也可以表示为幂级数：

However, exponentiation can also be expressed as a power series:

这个定义适用于实数、虚数，并且通过简单的推广也适用于复数。将此定义应用于纯虚数iϕ，我们得到：

This definition works with real numbers, with imaginary numbers and, by a simple extension, with complex numbers. Applying this definition to a pure imaginary number, iϕ, we get:

通过重新排列项，我们可以证明这等价于：

By rearranging terms, we can show that this is equivalent to:

您可以在http://en.wikipedia.org/wiki/Euler's_formula查看推导过程。

You can see the derivation at http://en.wikipedia.org/wiki/Euler’s_formula.

该公式表明e ^iϕ是一个模为 1 的复数；如果将其视为复平面上的一个点，它始终位于单位圆上。如果将其视为一个向量，则该向量与 x 轴正方向之间的夹角（以弧度为单位）即为辐角ϕ。

This formula implies that e^iϕ is a complex number with magnitude 1; if you think of it as a point in the complex plane, it is always on the unit circle. And if you think of it as a vector, the angle in radians between the vector and the positive x-axis is the argument ϕ.

当指数为复数时，我们有：

In the case where the exponent is a complex number, we have:

其中A是表示振幅的实数，e ^iϕ是表示角度的单位复数。

where A is a real number that indicates amplitude and e^iϕ is a unit complex number that indicates angle.

expNumPy 提供了一个可以处理复数的版本：

NumPy provides a version of exp that works with complex numbers:

>>> phi = 1.5 
>>> z = np.exp(1j * phi) 
>>> z 
(0.0707+0.997j)

>>> phi = 1.5
>>> z = np.exp(1j * phi)
>>> z
(0.0707+0.997j)

Python 使用 0j来表示虚数单位，而不是i0。以 0 结尾的数字j被认为是虚数，所以1j0 就是i。

Python uses j to represent the imaginary unit, rather than i. A number ending in j is considered imaginary, so 1j is just i.

当参数为np.exp虚数或复数时，结果为复数；具体来说，是一个复数np.complex128，它由两个 64 位浮点数表示。在本例中，结果为0.0707+0.997j。

When the argument to np.exp is imaginary or complex, the result is a complex number; specifically, an np.complex128, which is represented by two 64-bit floating-point numbers. In this example, the result is 0.0707+0.997j.

复数具有real以下属性imag：

Complex numbers have attributes real and imag:

>>> z.real
0.0707
>>> z.imag 
0.997

>>> z.real
0.0707
>>> z.imag
0.997

要获取幅度，您可以使用内置函数abs或np.absolute：

To get the magnitude, you can use the built-in function abs or np.absolute:

>>> abs(z)
1.0
>>> np.absolute(z) 
1.0

>>> abs(z)
1.0
>>> np.absolute(z)
1.0

要获取角度，您可以使用np.angle：

To get the angle, you can use np.angle:

>>> np.angle(z) 
1.5

>>> np.angle(z)
1.5

这个例子证实了是一个复数，其大小为 1，角度为ϕ弧度。

This example confirms that is a complex number with magnitude 1 and angle ϕ radians.

复杂信号

Complex Signals

如果是时间的函数，那么 也是时间的函数。具体来说：

If is a function of time, is also a function of time. Specifically:

该函数描述的是一个随时间变化的量，因此它是一个信号。具体来说，它是一个复指数信号。

This function describes a quantity that varies in time, so it is a signal. Specifically, it is a complex exponential signal.

在信号频率恒定的特殊情况下，结果是复正弦波：

In the special case where the frequency of the signal is constant, is and the result is a complex sinusoid:

更一般地，信号可能从相位偏移ϕ ₀开始，从而产生：

Or more generally, the signal might start at a phase offset ϕ₀, yielding:

thinkdsp提供此信号的实现ComplexSinusoid：

thinkdsp provides an implementation of this signal, ComplexSinusoid:

class ComplexSinusoid(Sinusoid):
 
   def evaluate(self, ts):
        相位 = PI2 * 自我频率 * 时间 + 自我偏移
        ys = self.amp * np.exp(1j * phases)
        返回 ys

class ComplexSinusoid(Sinusoid):
 
   def evaluate(self, ts):
        phases = PI2 * self.freq * ts + self.offset
        ys = self.amp * np.exp(1j * phases)
        return ys

ComplexSinusoid继承__init__自。它提供了一个几乎与相同的Sinusoid版本；唯一的区别是它使用而不是。evaluateSinusoid.evaluatenp.expnp.sin

ComplexSinusoid inherits __init__ from Sinusoid. It provides a version of evaluate that is almost identical to Sinusoid.evaluate; the only difference is that it uses np.exp instead of np.sin.

结果是一个包含复数的 NumPy 数组：

The result is a NumPy array of complex numbers:

>>> signal = thinkdsp.ComplexSinusoid(freq=1, amp=0.6, offset=1) 
>>> wave = signal.make_wave(duration=1, framerate=4) 
>>> wave.ys 
[ 0.324+0.505j -0.505+0.324j -0.324-0.505j 0.505-0.324j]

>>> signal = thinkdsp.ComplexSinusoid(freq=1, amp=0.6, offset=1)
>>> wave = signal.make_wave(duration=1, framerate=4)
>>> wave.ys
[ 0.324+0.505j -0.505+0.324j -0.324-0.505j  0.505-0.324j]

该信号的频率为每秒 1 个周期；幅度为 0.6（单位未指定）；相位偏移为 1 弧度。

The frequency of this signal is 1 cycle per second; the amplitude is 0.6 (in unspecified units); and the phase offset is 1 radian.

此示例在 0 到 1 秒之间等距选取四个位置对信号进行评估。所得样本为复数。

This example evaluates the signal at four places equally spaced between 0 and 1 second. The resulting samples are complex numbers.

合成问题

The Synthesis Problem

就像我们处理真实正弦波一样，我们可以通过将不同频率的复正弦波相加来创建复合信号。这就引出了复合合成问题的复数版本：给定每个复分量的频率和幅度，我们如何评估信号？

Just as we did with real sinusoids, we can create compound signals by adding up complex sinusoids with different frequencies. And that brings us to the complex version of the synthesis problem: given the frequency and amplitude of each complex component, how do we evaluate the signal?

最简单的解决方法是创建ComplexSinusoid对象并将它们相加：

The simplest solution is to create ComplexSinusoid objects and add them up:

def synthesize1(amps, fs, ts):
    components = [thinkdsp.ComplexSinusoid(freq, amp)
                  对于 zip(amps, fs) 中的每个 amp 和 freq]
    signal = thinkdsp.SumSignal(*components)
    ys = signal.evaluate(ts)
    返回 ys

def synthesize1(amps, fs, ts):
    components = [thinkdsp.ComplexSinusoid(freq, amp)
                  for amp, freq in zip(amps, fs)]
    signal = thinkdsp.SumSignal(*components)
    ys = signal.evaluate(ts)
    return ys

这个函数几乎与“Synthesis”synthesize1中的函数完全相同；唯一的区别是我用替换了。CosSignalComplexSinusoid

This function is almost identical to synthesize1 in “Synthesis”; the only difference is that I replaced CosSignal with ComplexSinusoid.

举个例子：

Here’s an example:

amps = np.array([0.6, 0.25, 0.1, 0.05])
fs = [100, 200, 300, 400]
帧率 = 11025
ts = np.linspace(0, 1, framerate)
ys = synthesize1(amps, fs, ts)

amps = np.array([0.6, 0.25, 0.1, 0.05])
fs = [100, 200, 300, 400]
framerate = 11025
ts = np.linspace(0, 1, framerate)
ys = synthesize1(amps, fs, ts)

结果是：

The result is:

[ 1.000 +0.000e+00j 0.995 +9.093e-02j 0.979 +1.803e-01j ...,
  0.979 -1.803e-01j 0.995 -9.093e-02j 1.000 -5.081e-15j]

[ 1.000 +0.000e+00j  0.995 +9.093e-02j  0.979 +1.803e-01j ...,
  0.979 -1.803e-01j  0.995 -9.093e-02j  1.000 -5.081e-15j]

从最基本的层面来说，复杂信号是一串复数。但我们应该如何解释它呢？我们对实际信号有一些直觉：它们代表随时间变化的量；例如，声音信号代表气压的变化。但我们在现实世界中测量的任何事物都不会产生复数。

At the lowest level, a complex signal is a sequence of complex numbers. But how should we interpret it? We have some intuition for real signals: they represent quantities that vary in time; for example, a sound signal represents changes in air pressure. But nothing we measure in the world yields complex numbers.

那么，什么是复杂信号？我没有令人满意的答案。我能提供的最好答案是两个不太令人满意的答案：

So what is a complex signal? I don’t have a satisfying answer to this question. The best I can offer is two unsatisfying answers:

1. 复杂信号是一种数学抽象概念，对计算和分析很有用，但它并不直接对应于现实世界中的任何事物。

2. A complex signal is a mathematical abstraction that is useful for computation and analysis, but it does not correspond directly with anything in the real world.

3. 如果你愿意，你可以把复杂信号想象成一个复数序列，其中包含两个信号作为其实部和虚部。

4. If you like, you can think of a complex signal as a sequence of complex numbers that contains two signals as its real and imaginary parts.

从第二个角度来看，我们可以将前面的信号分解为实部和虚部：

Taking the second point of view, we can split the previous signal into its real and imaginary parts:

n = 500
thinkplot.plot(ts[:n], ys[:n].real, label='real')
thinkplot.plot(ts[:n], ys[:n].imag, label='imag')

n = 500
thinkplot.plot(ts[:n], ys[:n].real, label='real')
thinkplot.plot(ts[:n], ys[:n].imag, label='imag')

图 7-1显示了结果的一部分。实部是余弦波之和；虚部是正弦波之和。虽然波形看起来不同，但它们包含相同比例的相同频率分量。在我们听来，它们听起来是一样的（通常情况下，我们听不出相位偏移）。

Figure 7-1 shows a segment of the result. The real part is a sum of cosines; the imaginary part is a sum of sines. Although the waveforms look different, they contain the same frequency components in the same proportions. To our ears, they sound the same (in general, we don’t hear phase offsets).

图 7-1.复杂正弦波混合物的实部和虚部。

基质合成

Synthesis with Matrices

正如我们在“使用数组进行综合”一节中所看到的，我们也可以用矩阵乘法来表示综合问题：

As we saw in “Synthesis with Arrays”, we can also express the synthesis problem in terms of matrix multiplication:

PI2 = np.pi * 2

def synthesize2(amps, fs, ts):
    args = np.outer(ts, fs)
    M = np.exp(1j * PI2 * args)
    ys = np.dot(M, amps)
    返回 ys

PI2 = np.pi * 2

def synthesize2(amps, fs, ts):
    args = np.outer(ts, fs)
    M = np.exp(1j * PI2 * args)
    ys = np.dot(M, amps)
    return ys

再次强调，amps这是一个包含一系列振幅的 NumPy 数组。

Again, amps is a NumPy array that contains a sequence of amplitudes.

fs是一个包含各分量频率的序列。ts包含我们将评估信号的时间点。

fs is a sequence containing the frequencies of the components. ts contains the times where we will evaluate the signal.

argsts包含和 的外积fs，其中 沿ts行向下，fs沿列向上（您可能需要回顾图 6-1）。

args contains the outer product of ts and fs, with the ts running down the rows and the fs running across the columns (you might want to refer back to Figure 6-1).

矩阵的每一列M都包含一个具有特定频率的复正弦波，该复正弦波在一系列处进行评估ts。

Each column of matrix M contains a complex sinusoid with a particular frequency, evaluated at a sequence of ts.

当我们乘以M振幅时，结果是一个向量，其元素对应于ts；每个元素都是在特定时间计算的几个复正弦波之和。

When we multiply M by the amplitudes, the result is a vector whose elements correspond to the ts; each element is the sum of several complex sinusoids, evaluated at a particular time.

以下是上一节中的例子：

Here’s the example from the previous section again:

>>> ys = synthesize2(amps, fs, ts) 
>>> ys
[ 1.000 +0.000e+00j 0.995 +9.093e-02j 0.979 +1.803e-01j ...,
  0.979 -1.803e-01j 0.995 -9.093e-02j 1.000 -5.081e-15j]

>>> ys = synthesize2(amps, fs, ts)
>>> ys
[ 1.000 +0.000e+00j  0.995 +9.093e-02j  0.979 +1.803e-01j ...,
  0.979 -1.803e-01j  0.995 -9.093e-02j  1.000 -5.081e-15j]

结果是一样的。

The result is the same.

在这个例子中，振幅是实数，但也可以是复数。复振幅会对结果产生什么影响？记住，我们可以用两种方式理解复数：实部和虚部之和，或者实振幅与复指数的乘积。使用第二种解释，我们可以看到当复振幅乘以复正弦波时会发生什么。对于每个频率f，我们有：

In this example the amplitudes are real, but they could also be complex. What effect does a complex amplitude have on the result? Remember that we can think of a complex number in two ways: either the sum of a real and an imaginary part, , or the product of a real amplitude and a complex exponential, . Using the second interpretation, we can see what happens when we multiply a complex amplitude by a complex sinusoid. For each frequency, f, we have:

乘以A将振幅乘以A，并加上相位偏移ϕ ₀。

Multiplying by multiplies the amplitude by A and adds the phase offset ϕ₀.

我们可以通过运行前面的例子，并使用复振幅来验证这个说法：

We can test that claim by running the previous example with complex amplitudes:

φ = 1.5
amps2 = amps * np.exp(1j * phi)
ys2 = synthesize2(amps2, fs, ts)

thinkplot.plot(ts[:n], ys.real[:n])
thinkplot.plot(ts[:n], ys2.real[:n])

phi = 1.5
amps2 = amps * np.exp(1j * phi)
ys2 = synthesize2(amps2, fs, ts)

thinkplot.plot(ts[:n], ys.real[:n])
thinkplot.plot(ts[:n], ys2.real[:n])

由于amps是一个实数数组，乘以 会np.exp(1j * phi)产生一个复数数组，其相位偏移phi弧度为 ，幅度与 相同amps。

Since amps is an array of reals, multiplying by np.exp(1j * phi) yields an array of complex numbers with phase offset phi radians, and the same magnitudes as amps.

图 7-2显示了具有不同相位偏移的波形。每个频率分量都偏移了大约四分之一周期。但不同频率的分量具有不同的周期；因此，每个分量在时间上的偏移量也不同。当我们把这些分量叠加起来时，得到的波形看起来就不一样了。

Figure 7-2 shows waveforms with different phase offsets. With each frequency component gets shifted by about a quarter of a cycle. But components with different frequencies have different cycles; as a result, each component is shifted by a different amount in time. When we add up the components, the resulting waveforms look different.

现在我们有了更通用的合成问题解决方案——一个可以处理复杂振幅的解决方案——我们就可以着手解决分析问题了。

Now that we have the more general solution to the synthesis problem—one that handles complex amplitudes—we are ready for the analysis problem.

图 7-2.两个相差一个相位偏移的复信号的实部。

分析问题

The Analysis Problem

分析问题是合成问题的逆问题：给定一个样本序列y，并且知道构成信号的频率，我们能否计算出分量的复振幅a？

The analysis problem is the inverse of the synthesis problem: given a sequence of samples, y, and knowing the frequencies that make up the signal, can we compute the complex amplitudes of the components, a?

正如我们在“分析”部分所看到的，我们可以通过构造综合矩阵M并求解线性方程组 来解决这个问题，其中a为参数：

As we saw in “Analysis”, we can solve this problem by forming the synthesis matrix, M, and solving the system of linear equations, , for a:

def analyze1(ys, fs, ts):
    args = np.outer(ts, fs)
    M = np.exp(1j * PI2 * args)
    amps = np.linalg.solve(M, ys)
    回流放大器

def analyze1(ys, fs, ts):
    args = np.outer(ts, fs)
    M = np.exp(1j * PI2 * args)
    amps = np.linalg.solve(M, ys)
    return amps

analyze1它接受一个（可能复杂的）波数组ys，一个实频率序列，fs以及一个实时间序列ts。它返回一个复振幅序列amps。

analyze1 takes a (possibly complex) wave array, ys, a sequence of real frequencies, fs, and a sequence of real times, ts. It returns a sequence of complex amplitudes, amps.

继续之前的例子，我们可以确认它analyze1恢复了我们最初的振幅。为了使线性系统求解器正常工作，M必须是正方形，所以我们需要ys、fs和ts具有相同的长度。我将通过将ys和ts切片到 的长度来确保这一点fs：

Continuing the previous example, we can confirm that analyze1 recovers the amplitudes we started with. For the linear system solver to work, M has to be square, so we need ys, fs and ts, to have the same length. I’ll ensure that by slicing ys and ts down to the length of fs:

>>> n = len(fs) 
>>> amps2 = analyze1(ys[:n], fs, ts[:n]) 
>>> amps2 
[ 0.60+0.j 0.25-0.j 0.10+0.j 0.05-0.j]

>>> n = len(fs)
>>> amps2 = analyze1(ys[:n], fs, ts[:n])
>>> amps2
[ 0.60+0.j  0.25-0.j  0.10+0.j  0.05-0.j]

这些与我们最初得到的振幅大致相同，尽管由于浮点误差，每个分量都有一个很小的虚部。

These are approximately the amplitudes we started with, although each component has a small imaginary part due to floating-point errors.

高效分析

Efficient Analysis

遗憾的是，求解线性方程组速度很慢。对于离散余弦变换（DCT），我们通过选择合适的矩阵使得矩阵正交来加快计算速度。fs这样，矩阵的逆矩阵就是矩阵的转置矩阵，我们就可以通过矩阵乘法来计算DCT及其逆矩阵。tsMMM

Unfortunately, solving a linear system of equations is slow. For the DCT, we were able to speed things up by choosing fs and ts so that M is orthogonal. That way, the inverse of M is the transpose of M, and we can compute both the DCT and inverse DCT by matrix multiplication.

我们对 DFT 也做同样的处理，但有一个小小的改动。由于M是复数，我们需要它是酉矩阵而不是正交矩阵，这意味着 的逆矩阵M是 的共轭转置矩阵M，我们可以通过转置矩阵并对每个元素的虚部取反来计算共轭转置矩阵。参见http://en.wikipedia.org/wiki/Unitary_matrix。

We’ll do the same thing for the DFT, with one small change. Since M is complex, we need it to be unitary, rather than orthogonal, which means that the inverse of M is the conjugate transpose of M, which we can compute by transposing the matrix and negating the imaginary part of each element. See http://en.wikipedia.org/wiki/Unitary_matrix.

这些 NumPy 方法conj可以实现我们想要的功能。以下是计算分量的transpose代码：M

The NumPy methods conj and transpose do what we want. Here’s the code that computes M for components:

N = 4
ts = np.arange(N) / N
fs = np.arange(N)
args = np.outer(ts, fs)
M = np.exp(1j * PI2 * args)

N = 4
ts = np.arange(N) / N
fs = np.arange(N)
args = np.outer(ts, fs)
M = np.exp(1j * PI2 * args)

如果M是酉矩阵，则，其中是M的共轭转置，I是单位矩阵。我们可以这样检验M是否为酉矩阵：

If M is unitary, , where is the conjugate transpose of M, and I is the identity matrix. We can test whether M is unitary like this:

MstarM = M.conj().transpose().dot(M)

MstarM = M.conj().transpose().dot(M)

在浮点误差容限内，结果是，因此M除了一个额外的因子N之外是单位矩阵，类似于我们在 DCT 中发现的额外因子 2。

The result, within the tolerance of floating-point error, is , so M is unitary except for an extra factor of N, similar to the extra factor of 2 we found with the DCT.

我们可以利用这个结果编写一个速度更快的版本analyze1：

We can use this result to write a faster version of analyze1:

def analyze2(ys, fs, ts):
    args = np.outer(ts, fs)
    M = np.exp(1j * PI2 * args)
    amps = M.conj().transpose().dot(ys) / N
    回流放大器

def analyze2(ys, fs, ts):
    args = np.outer(ts, fs)
    M = np.exp(1j * PI2 * args)
    amps = M.conj().transpose().dot(ys) / N
    return amps

fs并使用适当的和值进行测试ts：

And test it with appropriate values of fs and ts:

N = 4
amps = np.array([0.6, 0.25, 0.1, 0.05])
fs = np.arange(N)
ts = np.arange(N) / N
ys = synthesize2(amps, fs, ts)
amps3 = analyze2(ys, fs, ts)

N = 4
amps = np.array([0.6, 0.25, 0.1, 0.05])
fs = np.arange(N)
ts = np.arange(N) / N
ys = synthesize2(amps, fs, ts)
amps3 = analyze2(ys, fs, ts)

同样，结果在浮点运算的误差范围内是正确的：

Again, the result is correct within the tolerance of floating-point arithmetic:

[ 0.60+0.j 0.25+0.j 0.10-0.j 0.05-0.j]

[ 0.60+0.j  0.25+0.j  0.10-0.j  0.05-0.j]

密度泛函理论

DFT

作为函数，analyze2它很难使用，因为它只有在正确选择 `x`fs和 `y`ts参数时才能正常工作。因此，我将重写它，使其只接受 `x`ys和 `y` 参数，并计算freq它们ts自身。

As a function, analyze2 would be hard to use because it only works if fs and ts are chosen correctly. Instead, I will rewrite it to take just ys and compute freq and ts itself.

首先，我将编写一个函数来计算综合矩阵M：

First, I’ll make a function to compute the synthesis matrix, M:

def synthesis_matrix(N):
    ts = np.arange(N) / N
    fs = np.arange(N)
    args = np.outer(ts, fs)
    M = np.exp(1j * PI2 * args)
    返回 M

def synthesis_matrix(N):
    ts = np.arange(N) / N
    fs = np.arange(N)
    args = np.outer(ts, fs)
    M = np.exp(1j * PI2 * args)
    return M

ys然后我将编写一个接受参数并返回结果的函数amps：

Then I’ll write the function that takes ys and returns amps:

def analyze3(ys):
    N = len(ys)
    M = synthesis_matrix(N)
    amps = M.conj().transpose().dot(ys) / N
    回流放大器

def analyze3(ys):
    N = len(ys)
    M = synthesis_matrix(N)
    amps = M.conj().transpose().dot(ys) / N
    return amps

我们快完成了；analyze3计算结果与DFT非常接近，只有一个区别。DFT的传统定义不进行除法运算N：

We are almost done; analyze3 computes something very close to the DFT, with one difference. The conventional definition of the DFT does not divide by N:

def dft(ys):
    N = len(ys)
    M = synthesis_matrix(N)
    amps = M.conj().transpose().dot(ys)
    回流放大器

def dft(ys):
    N = len(ys)
    M = synthesis_matrix(N)
    amps = M.conj().transpose().dot(ys)
    return amps

现在我们可以确认，我的版本与以下版本的结果相同np.fft.fft：

Now we can confirm that my version yields the same result as np.fft.fft:

>>> dft(ys) 
[ 2.4+0.j 1.0+0.j 0.4-0.j 0.2-0.j]

>>> dft(ys)
[ 2.4+0.j  1.0+0.j  0.4-0.j  0.2-0.j]

结果接近于amps * N。以下是版本np.fft：

The result is close to amps * N. And here’s the version in np.fft:

>>> np.fft.fft(ys) 
[ 2.4+0.j 1.0+0.j 0.4-0.j 0.2-0.j]

>>> np.fft.fft(ys)
[ 2.4+0.j  1.0+0.j  0.4-0.j  0.2-0.j]

在浮点误差范围内，它们是相同的。

They are the same, within floating-point error.

逆DFT几乎相同，只是我们不需要转置和共轭M，现在我们需要除以N：

The inverse DFT is almost the same, except we don’t have to transpose and conjugate M, and now we have to divide through by N:

def idft(ys):
    N = len(ys)
    M = synthesis_matrix(N)
    安培数 = M.dot(ys) / N
    回流放大器

def idft(ys):
    N = len(ys)
    M = synthesis_matrix(N)
    amps = M.dot(ys) / N
    return amps

最后，我们可以确认dft(idft(amps))产量为amps：

Finally, we can confirm that dft(idft(amps)) yields amps:

>>> ys = idft(amps) 
>>> dft(ys) 
[ 0.60+0.j 0.25+0.j 0.10-0.j 0.05-0.j]

>>> ys = idft(amps)
>>> dft(ys)
[ 0.60+0.j  0.25+0.j  0.10-0.j  0.05-0.j]

如果可以回到过去，我可能会修改DFT的定义，使其除以N，而逆DFT则不除以N。这样会更符合我对合成和分析问题的阐述。

If I could go back in time, I might change the definition of the DFT so it divides by N and the inverse DFT doesn’t. That would be more consistent with my presentation of the synthesis and analysis problems.

或者我可以修改定义，使两个运算都除以。这样，DFT 和逆 DFT 就会更加对称。

Or I might change the definition so that both operations divide through by . Then the DFT and inverse DFT would be more symmetric.

但我现在还回不去过去，所以我们只能接受这种略显怪异的惯例。不过实际上，这倒无关紧要。

But I can’t go back in time (yet!), so we’re stuck with a slightly weird convention. For practical purposes it doesn’t really matter.

密度泛函理论是周期性的。

The DFT Is Periodic

本章中，我以矩阵乘法的形式介绍了DFT。我们计算合成矩阵M和分析矩阵。当乘以波阵列y时，结果的每个元素都是M和y的某一行的乘积，我们可以将其写成求和的形式：

In this chapter I presented the DFT in the form of matrix multiplication. We compute the synthesis matrix, M, and the analysis matrix, . When we multiply by the wave array, y, each element of the result is the product of a row from and y, which we can write in the form of a summation:

其中k是频率的索引，范围从 0 到，n是时间的索引，范围从 0 到。因此，是y的 DFT 的第k个元素。

where k is an index of frequency from 0 to and n is an index of time from 0 to . So is the kth element of the DFT of y.

通常情况下，我们会对N 个k值（从 0 到n）计算这个求和式。我们也可以计算其他k值，但这没有意义，因为结果会重复。也就是说，k处的值与 k = 1或k = 2处的值相同，等等。

Normally we evaluate this summation for N values of k, running from 0 to . We could evaluate it for other values of k, but there is no point, because they start to repeat. That is, the value at k is the same as the value at or or , etc.

我们可以通过代入求和式来从数学上看出这一点：

We can see that mathematically by plugging into the summation:

由于指数部分包含求和，我们可以将其拆分为两部分：

Since there is a sum in the exponent, we can break it into two parts:

在第二项中，指数始终是 的整数倍，因此结果始终为 1，我们可以将其省略：

In the second term, the exponent is always an integer multiple of , so the result is always 1, and we can drop it:

我们可以看到，这个求和等价于。因此，DFT 是周期性的，周期为N。本章末尾的练习题要求你实现快速傅里叶变换 (FFT)，你会用到这个结果。

And we can see that this summation is equivalent to . So the DFT is periodic, with period N. You will need this result for one of the exercises at the end of this chapter, which asks you to implement the Fast Fourier Transform (FFT).

顺便一提，将DFT写成求和的形式有助于理解它的工作原理。如果你回顾一下“使用数组进行合成”中的图表，你会发现合成矩阵的每一列都是一个信号在一系列时间点上的值。分析矩阵是合成矩阵的（共轭）转置，因此每一行都是一个信号在一系列时间点上的值。

As an aside, writing the DFT in the form of a summation provides an insight into how it works. If you review the diagram in “Synthesis with Arrays”, you’ll see that each column of the synthesis matrix is a signal evaluated at a sequence of times. The analysis matrix is the (conjugate) transpose of the synthesis matrix, so each row is a signal evaluated at a sequence of times.

因此，每个求和项都是y与阵列中某个信号的相关性（参见“相关性作为点积” ）。也就是说，DFT 的每个元素都是一个相关性，它量化了波阵列y与特定频率下的复指数函数的相似性。

Therefore, each summation is the correlation of y with one of the signals in the array (see “Correlation as Dot Product”). That is, each element of the DFT is a correlation that quantifies the similarity of the wave array, y, and a complex exponential at a particular frequency.

实信号的离散傅里叶变换

DFT of Real Signals

Spectrum本章中的类thinkdsp基于计算“实数DFT”的函数np.ftt.rfft，也就是说，它处理的是实数信号。但本章介绍的DFT更具通用性，它适用于复数信号。

The Spectrum class in thinkdsp is based on np.ftt.rfft, which computes the “real DFT”; that is, it works with real signals. But the DFT as presented in this chapter is more general than that; it works with complex signals.

那么，当我们对实际信号应用“完整DFT”时会发生什么呢？让我们来看一个例子：

So what happens when we apply the “full DFT” to a real signal? Let’s look at an example:

signal = thinkdsp.SawtoothSignal(freq=500)
wave = signal.make_wave(duration=0.1, framerate=10000)
hs = dft(wave.ys)
amps = np.absolute(hs)

signal = thinkdsp.SawtoothSignal(freq=500)
wave = signal.make_wave(duration=0.1, framerate=10000)
hs = dft(wave.ys)
amps = np.absolute(hs)

这段代码生成一个频率为 500 Hz、采样率为 10 kHz 的锯齿波。hs它包含了该波的复数离散傅里叶变换 (DFT)；amps包含了每个频率处的振幅。但是这些振幅分别对应于哪个频率呢？如果我们观察函数体dft，会发现：

This code makes a sawtooth wave with frequency 500 Hz, sampled at frame rate 10 kHz. hs contains the complex DFT of the wave; amps contains the amplitude at each frequency. But what frequency do these amplitudes correspond to? If we look at the body of dft, we see:

fs = np.arange(N)

fs = np.arange(N)

人们很容易认为这些值就是正确的频率。问题在于它dft不知道采样率。DFT 假设波形的持续时间为 1 个时间单位，因此它认为采样率为每时间单位N 次。为了解释这些频率，我们必须将这些任意的时间单位转换回秒，如下所示：

It’s tempting to think that these values are the right frequencies. The problem is that dft doesn’t know the sampling rate. The DFT assumes that the duration of the wave is 1 time unit, so it thinks the sampling rate is N per time unit. In order to interpret the frequencies, we have to convert from these arbitrary time units back to seconds, like this:

fs = np.arange(N) * 帧率​​ / N

fs = np.arange(N) * framerate / N

经过这一更改，频率范围从 0 到实际帧速率 10 kHz。现在我们可以绘制频谱图了：

With this change, the range of frequencies is from 0 to the actual frame rate, 10 kHz. Now we can plot the spectrum:

thinkplot.plot(fs, amps)
thinkplot.config(xlabel='频率 (Hz)',
                 ylabel='振幅')

thinkplot.plot(fs, amps)
thinkplot.config(xlabel='frequency (Hz)', 
                 ylabel='amplitude')

图 7-3显示了 0 到 10 kHz 范围内每个频率分量的信号幅度。图的左半部分符合我们的预期：主频率为 500 Hz，谐波衰减如下。

Figure 7-3 shows the amplitude of the signal for each frequency component from 0 to 10 kHz. The left half of the figure is what we should expect: the dominant frequency is at 500 Hz, with harmonics dropping off like .

图 7-3.采样率为 10 kHz 的 500 Hz 锯齿波信号的 DFT。

但图中右半部分却令人惊讶。超过 5000 赫兹后，谐波的振幅再次开始增大，并在 9500 赫兹处达到峰值。这是怎么回事？

But the right half of the figure is a surprise. Past 5000 Hz, the amplitude of the harmonics starts growing again, peaking at 9500 Hz. What’s going on?

答案是：混叠。请记住，帧速率为 10,000 Hz 时，折叠频率为 5000 Hz。正如我们在“混叠”一节中看到的，5500 Hz 的分量与 4500 Hz 的分量无法区分。当我们计算 5500 Hz 处的 DFT 时，会得到与 4500 Hz 处相同的值。类似地，6000 Hz 处的值与 4000 Hz 处的值相同，依此类推。

The answer: aliasing. Remember that with frame rate 10,000 Hz, the folding frequency is 5000 Hz. As we saw in “Aliasing”, a component at 5500 Hz is indistinguishable from a component at 4500 Hz. When we evaluate the DFT at 5500 Hz, we get the same value as at 4500 Hz. Similarly, the value at 6000 Hz is the same as the one at 4000 Hz, and so on.

实际信号的离散傅里叶变换 (DFT) 关于折叠频率对称。由于超过折叠频率后没有额外信息，我们可以只计算 DFT 的前半部分，从而节省时间，而这正是我们np.fft.rfft所做的。

The DFT of a real signal is symmetric around the folding frequency. Since there is no additional information past this point, we can save time by evaluating only the first half of the DFT, and that’s exactly what np.fft.rfft does.

平滑

Smoothing

平滑处理是一种旨在消除信号中短期波动以揭示长期趋势的操作。例如，如果您绘制股票价格的每日变化图，图表会显得杂乱无章；平滑处理可以使价格随时间推移的总体趋势（上涨或下跌）更加清晰。

Smoothing is an operation that tries to remove short-term variations from a signal in order to reveal long-term trends. For example, if you were to plot daily changes in the price of a stock, it would look noisy; a smoothing operator might make it easier to see whether the price was generally going up or down over time.

常见的平滑算法是移动平均，它计算前n 个值的平均值，其中n为某个值。

A common smoothing algorithm is a moving average, which computes the mean of the previous n values, for some value of n.

例如，图 8-1显示了 Facebook 股票从 2012 年 5 月 17 日至 2015 年 12 月 8 日的每日收盘价。灰色线代表原始数据，深色线代表 30 天移动平均线。平滑处理可以消除极端波动，使长期趋势更容易显现。

For example, Figure 8-1 shows the daily closing price of Facebook stock from May 17, 2012 to December 8, 2015. The gray line is the raw data, and the darker line shows the 30-day moving average. Smoothing removes the most extreme changes and makes it easier to see long-term trends.

图 8-1. Facebook 股票的每日收盘价和 30 天移动平均线。

平滑操作也适用于音频信号。例如，我将从一个 440 Hz 的方波开始。正如我们在“方波”一节中看到的，方波的谐波衰减缓慢，因此它包含许多高频分量。

Smoothing operations also apply to sound signals. As an example, I’ll start with a square wave at 440 Hz. As we saw in “Square Waves”, the harmonics of a square wave drop off slowly, so it contains many high-frequency components.

首先，我将构建信号和两个波形：

First I’ll construct the signal and two waves:

signal = thinkdsp.SquareSignal(freq=440)
wave = signal.make_wave(duration=1, framerate=44100)
segment = wave.segment(duration=0.01)

signal = thinkdsp.SquareSignal(freq=440)
wave = signal.make_wave(duration=1, framerate=44100)
segment = wave.segment(duration=0.01)

wave是信号的一秒钟切片；segment是我将用于绘图的较短切片。

wave is a one-second slice of the signal; segment is a shorter slice I’ll use for plotting.

为了计算该信号的移动平均值，我将使用类似于“加窗”一节中介绍的窗口。之前我们使用汉明窗来避免信号开头和结尾处的不连续性导致的频谱泄漏。更一般地，我们可以使用窗口来计算波形中样本的加权和。

To compute the moving average of this signal, I’ll use a window similar to the ones in “Windowing”. Previously we used a Hamming window to avoid spectral leakage caused by discontinuity at the beginning and end of a signal. More generally, we can use windows to compute the weighted sum of samples in a wave.

例如，要计算移动平均值，我将创建一个包含 11 个元素的窗口，并对其进行归一化，使这些元素的总和为 1：

For example, to compute a moving average, I’ll create a window with 11 elements and normalize it so the elements add up to 1:

window = np.ones(11)
window /= sum(window)

window = np.ones(11)
window /= sum(window)

现在我可以通过将窗口大小乘以波形数组来计算前 11 个元素的平均值：

Now I can compute the average of the first 11 elements by multiplying the window by the wave array:

ys = segment.ys
N = len(ys)
padded = thinkdsp.zero_pad(window, N)
prod = padded * ys
sum(prod)

ys = segment.ys
N = len(ys)
padded = thinkdsp.zero_pad(window, N)
prod = padded * ys
sum(prod)

padded是在窗口末尾添加零使其长度与 相同的版本segment.ys。像这样添加零称为填充。

padded is a version of the window with zeros added to the end so it is the same length as segment.ys. Adding zeros like this is called padding.

prod是窗口和波形阵列的乘积。各元素乘积之和等于阵列前 11 个元素的平均值。由于这些元素均为 -1，因此它们的平均值为 -1。

prod is the product of the window and the wave array. The sum of the elementwise products is the average of the first 11 elements of the array. Since these elements are all –1, their average is –1.

为了计算移动平均值的下一个元素，我们滚动窗口，将 1 向右移动，并将末尾的一个 0 循环到开头。

To compute the next element of the moving average, we roll the window, which shifts the ones to the right and wraps one of the zeros from the end around to the beginning.

当我们把滚动窗口和波阵面相乘时，我们得到波阵面接下来 11 个元素的平均值，从第二个元素开始：

When we multiply the rolled window and the wave array, we get the average of the next 11 elements of the wave array, starting with the second:

rolling = np.roll(rolled, 1)
产品 = 滚动 * ys
sum(prod)

rolled = np.roll(rolled, 1)
prod = rolled * ys
sum(prod)

结果仍然是-1。

The result is –1 again.

我们可以用同样的方法计算其余元素。下面的函数将我们目前看到的代码封装在一个循环中，并将结果存储在一个数组中：

We can compute the rest of the elements the same way. The following function wraps the code we have seen so far in a loop and stores the results in an array:

def smooth(ys, window):
    N = len(ys)
    smoothed = np.zeros(N)
    padded = thinkdsp.zero_pad(window, N)
    卷边 = 填充

    for i in range(N):
        smoothed[i] = sum(rolled * ys)
        rolling = np.roll(rolled, 1)
    返回平滑

def smooth(ys, window):
    N = len(ys)
    smoothed = np.zeros(N)
    padded = thinkdsp.zero_pad(window, N)
    rolled = padded

    for i in range(N):
        smoothed[i] = sum(rolled * ys)
        rolled = np.roll(rolled, 1)
    return smoothed

smoothed是包含结果的数组；padded是包含窗口和足够多的零以达到长度的数组N；rolled是的副本，padded每次循环时都会向右移动一个元素。

smoothed is the array that will contain the results; padded is an array that contains the window and enough zeros to have length N; and rolled is a copy of padded that gets shifted to the right by one element each time through the loop.

在循环内部，我们ys乘以rolled来选择 11 个元素并将它们相加。

Inside the loop, we multiply ys by rolled to select 11 elements and add them up.

图 8-2显示了方波的处理结果。灰色线代表原始信号；深色线代表平滑后的信号。平滑后的信号在窗口前沿到达第一个过渡点时开始上升，并在窗口越过过渡点时趋于平缓。因此，过渡点不再那么突兀，拐点也更加柔和。聆听平滑后的信号，会感觉嗡嗡声减少，声音略微变得沉闷。

Figure 8-2 shows the result for a square wave. The gray line is the original signal; the darker line is the smoothed signal. The smoothed signal starts to ramp up when the leading edge of the window reaches the first transition, and levels off when the window crosses the transition. As a result, the transitions are less abrupt, and the corners less sharp. If you listen to the smoothed signal, it sounds less buzzy and slightly muffled.

图 8-2. 400 Hz 的方波信号（灰色）和 11 元素移动平均值。

卷积

Convolution

我们刚才执行的操作——将窗口函数应用于波形的每个重叠部分——称为卷积。

The operation we just performed—applying a window function to each overlapping segment of a wave—is called convolution.

卷积运算非常常见，以至于 NumPy 提供了一个比我的版本更简单、更快速的实现：

Convolution is such a common operation that NumPy provides an implementation that is simpler and faster than my version:

convolve = np.convolve(ys, window, mode='valid')
smooth2 = thinkdsp.Wave(convolved, framerate=wave.framerate)

convolved = np.convolve(ys, window, mode='valid')
smooth2 = thinkdsp.Wave(convolved, framerate=wave.framerate)

np.convolve计算波阵列和窗口的卷积。模式标志指示它仅在窗口和波阵列完全重叠时才计算值，因此当窗口的右边缘到达波阵列的末端时，计算就会停止。除此之外，结果与图 8-2valid相同。

np.convolve computes the convolution of the wave array and the window. The mode flag valid indicates that it should only compute values when the window and the wave array overlap completely, so it stops when the right edge of the window reaches the end of the wave array. Other than that, the result is the same as in Figure 8-2.

实际上，还有另一点不同。上一节中的循环实际上计算的是互相关：

Actually, there is one other difference. The loop in the previous section actually computes cross-correlation:

其中f是长度为N 的波阵列，g是窗口，✩ 是互相关的符号。为了计算结果的第 n个元素，我们将g向右移动，这就是索引为 的原因。

where f is a wave array with length N, g is the window, and ✩ is the symbol for cross-correlation. To compute the nth element of the result, we shift g to the right, which is why the index is .

卷积的定义略有不同：

The definition of convolution is slightly different:

符号*表示卷积。区别在于g的索引：m被取反了，因此求和运算会反向遍历g的元素（假设负索引会循环到数组的末尾）。

The symbol * represents convolution. The difference is in the index of g: m has been negated, so the summation iterates the elements of g backward (assuming that negative indices wrap around to the end of the array).

由于本例中使用的窗口是对称的，因此互相关和卷积运算会得到相同的结果。使用其他窗口时，我们需要更加谨慎。

Because the window we used in this example is symmetric, cross-correlation and convolution yield the same result. When we use other windows, we will have to be more careful.

你可能会好奇为什么卷积的定义是这样，窗口的应用方式似乎与预期相反。原因有二：

You might wonder why convolution is defined like this, with the window applied in a way that seems backward. There are two reasons:

• 这个定义自然而然地出现在几个应用中，特别是信号处理系统的分析，这是第 10 章的主题。

• This definition comes up naturally for several applications, especially analysis of signal-processing systems, which is the topic of Chapter 10.

• 此外，这个定义是卷积定理的基础，我们很快就会讲到它。

• Also, this definition is the basis of the Convolution Theorem, coming up very soon.

最后，给那些了解太多的人提个醒：到目前为止，我还没有区分卷积和循环卷积。我们稍后会讲到。

Finally, a note for people who know too much: in the presentation so far I have not distinguished between convolution and circular convolution. We’ll get to it.

频域

The Frequency Domain

平滑处理使方波信号的过渡更加平缓，声音也略微变得沉闷。我们来看看这种操作对频谱的影响。首先，我将绘制原始波形的频谱：

Smoothing makes the transitions in a square signal less abrupt, and makes the sound slightly muffled. Let’s see what effect this operation has on the spectrum. First I’ll plot the spectrum of the original wave:

spectrum = wave.make_spectrum()
spectrum.plot(color=GRAY)

spectrum = wave.make_spectrum()
spectrum.plot(color=GRAY)

然后是平滑后的波形：

Then the smoothed wave:

convolve = np.convolve(wave.ys, window, mode='same')
smooth = thinkdsp.Wave(convolved, framerate=wave.framerate)
spectrum2 = smooth.make_spectrum()
spectrum2.plot()

convolved = np.convolve(wave.ys, window, mode='same')
smooth = thinkdsp.Wave(convolved, framerate=wave.framerate)
spectrum2 = smooth.make_spectrum()
spectrum2.plot()

模式标志same表明结果的长度应与输入的长度相同。在本例中，结果会包含一些“溢出”的值，但目前这没关系。

The mode flag same indicates that the result should have the same length as the input. In this example, it will include a few values that “wrap around”, but that’s OK for now.

图 8-3显示了结果。基频几乎不变；前几个谐波被衰减，高次谐波几乎完全消失。因此，平滑处理的效果类似于低通滤波器，我们在“频谱”和“粉红噪声”部分已经讨论过。

Figure 8-3 shows the result. The fundamental frequency is almost unchanged; the first few harmonics are attenuated, and the higher harmonics are almost eliminated. So smoothing has the effect of a low-pass filter, which we saw in “Spectrums” and “Pink Noise”.

图 8-3.平滑前后方波的频谱。

为了了解每个分量衰减了多少，我们可以计算两个光谱的比值：

To see how much each component has been attenuated, we can compute the ratio of the two spectrums:

安培数 = 频谱.安培数
amps2 = spectrum2.amps
比率 = 安培数² / 安培数    
ratio[amps<560] = 0
thinkplot.plot(ratio)

amps = spectrum.amps
amps2 = spectrum2.amps
ratio = amps2 / amps    
ratio[amps<560] = 0
thinkplot.plot(ratio)

ratio是平滑前后振幅的比值。当amps较小时，该比值可能很大且噪声较大，因此为了简化计算，除了谐波部分外，我将该比值设为 0。

ratio is the ratio of the amplitude before and after smoothing. When amps is small, this ratio can be big and noisy, so for simplicity I set the ratio to 0 except where the harmonics are.

图 8-4显示了结果。正如预期的那样，低频时比值很高，并在接近 4000 Hz 的截止频率处下降。但还有一个我们意想不到的特征：在截止频率以上，比值在 0 到 0.2 之间波动。这是怎么回事？

Figure 8-4 shows the result. As expected, the ratio is high for low frequencies and drops off at a cutoff frequency near 4000 Hz. But there is another feature we did not expect: above the cutoff, the ratio bounces around between 0 and 0.2. What’s up with that?

图 8-4.方波平滑前后的频谱比。

卷积定理

The Convolution Theorem

答案是卷积定理。用数学公式表述如下：

The answer is the Convolution Theorem. Stated mathematically:

其中f是波阵列，g是窗口。换句话说，卷积定理指出，如果我们对f和g进行卷积，然后计算它们的离散傅里叶变换 (DFT)，我们得到的结果与先分别计算f和g的 DFT ，然后将结果逐元素相乘的结果相同。

where f is a wave array and g is a window. In words, the Convolution Theorem says that if we convolve f and g, and then compute the DFT, we get the same answer as when computing the DFT of f and g, and then multiplying the results elementwise.

当我们对波形进行卷积等运算时，我们称是在时域中进行操作，因为波形是时间的函数。当我们对离散傅里叶变换（DFT）进行乘法等运算时，我们称是在频域中进行操作，因为DFT是频率的函数。

When we apply an operation like convolution to a wave, we say we are working in the time domain, because the wave is a function of time. When we apply an operation like multiplication to the DFT, we are working in the frequency domain, because the DFT is a function of frequency.

利用这些术语，我们可以更简洁地表述卷积定理：

Using these terms, we can state the Convolution Theorem more concisely:

时域中的卷积对应于频域中的乘法。

Convolution in the time domain corresponds to multiplication in the frequency domain.

这就解释了图 8-4 的含义，因为当我们对一个波和一个窗函数进行卷积时，实际上是将波的频谱与窗函数的频谱相乘。为了理解其工作原理，我们可以计算窗函数的 DFT：

And that explains Figure 8-4, because when we convolve a wave and a window, we multiply the spectrum of the wave with the spectrum of the window. To see how that works, we can compute the DFT of the window:

padded = zero_pad(window, N)
dft_window = np.fft.rfft(padded)
thinkplot.plot(abs(dft_window))

padded = zero_pad(window, N)
dft_window = np.fft.rfft(padded)
thinkplot.plot(abs(dft_window))

padded包含平滑窗口，用零填充使其长度与相同wave；dft_window包含的 DFT padded。

padded contains the smoothing window, padded with zeros to be the same length as wave; dft_window contains the DFT of padded.

图 8-5显示了结果，以及我们在上一节中计算的比值。这些比值恰好就是振幅dft_window。数学表达式为：

Figure 8-5 shows the result, along with the ratios we computed in the previous section. The ratios are exactly the amplitudes in dft_window. Mathematically:

在此上下文中，窗口的离散傅里叶变换（DFT）被称为滤波器。对于时域中的任何卷积窗口，在频域中都存在对应的滤波器。同样，对于频域中任何可以用元素级乘法表示的滤波器，也存在对应的窗口。

In this context, the DFT of a window is called a filter. For any convolution window in the time domain, there is a corresponding filter in the frequency domain. And for any filter that can be expressed by elementwise multiplication in the frequency domain, there is a corresponding window.

图 8-5.方波平滑前后的频谱比，以及平滑窗口的 DFT。

高斯滤波器

Gaussian Filter

我们在上一节中使用的移动平均窗口是一个低通滤波器，但它并不是一个很好的低通滤波器。离散傅里叶变换 (DFT) 的频谱起初下降很快，但随后会出现波动。这些波动被称为旁瓣，它们的出现是因为移动平均窗口类似于方波，因此其频谱包含高频谐波，而这些谐波的衰减速度相对较慢，与 成正比。

The moving average window we used in the previous section is a low-pass filter, but it is not a very good one. The DFT drops off steeply at first, but then it bounces around. Those bounces are called sidelobes, and they are there because the moving average window is like a square wave, so its spectrum contains high-frequency harmonics that drop off proportionally to , which is relatively slow.

使用高斯窗口可以做得更好。SciPy 提供了计算许多常用卷积窗口的函数，包括gaussian：

We can do better with a Gaussian window. SciPy provides functions that compute many common convolution windows, including gaussian:

高斯分布 = scipy.signal.gaussian(M=11, std=2)
高斯分布 /= ∑(高斯分布)

gaussian = scipy.signal.gaussian(M=11, std=2)
gaussian /= sum(gaussian)

M是窗口中元素的数量；std是用于计算该窗口的高斯分布的标准差。图 8-6显示了窗口的形状。它是高斯“钟形曲线”的离散近似。图中还显示了前一个例子中的移动平均窗口，它有时被称为“车厢窗口”，因为它看起来像一个矩形的火车车厢。

M is the number of elements in the window; std is the standard deviation of the Gaussian distribution used to compute it. Figure 8-6 shows the shape of the window. It is a discrete approximation of the Gaussian “bell curve”. The figure also shows the moving average window from the previous example, which is sometimes called a boxcar window because it looks like a rectangular railway car.

图 8-6.箱型窗和高斯窗。

我再次使用这个窗口运行了前面几节中的计算，并生成了图 8-7，该图显示了平滑前后光谱的比率以及高斯窗口的 DFT。

I ran the computations from the previous sections again with this window and generated Figure 8-7, which shows the ratio of the spectrums before and after smoothing, along with the DFT of the Gaussian window.

作为低通滤波器，高斯平滑比简单的移动平均效果更好。在比率下降后，它仍然保持在较低水平，几乎没有我们之前使用箱形窗口时看到的旁瓣。因此，它在滤除高频信号方面表现更佳。

As a low-pass filter, Gaussian smoothing is better than a simple moving average. After the ratio drops off, it stays low, with almost none of the sidelobes we saw with the boxcar window. So it does a better job of cutting off the higher frequencies.

它之所以表现如此出色，是因为高斯曲线的 DFT 也是一条高斯曲线。因此，该比率随 成比例下降，这比 下降得快得多。

The reason it does so well is that the DFT of a Gaussian curve is also a Gaussian curve. So the ratio drops off in proportion to , which is much faster than .

图 8-7.高斯平滑前后光谱的比值，以及窗口的 DFT。

高效卷积

Efficient Convolution

FFT 之所以如此重要，其中一个原因是它与卷积定理相结合，提供了一种计算卷积、互相关和自相关的有效方法。

One of the reasons the FFT is such an important algorithm is that, combined with the Convolution Theorem, it provides an efficient way to compute convolution, cross-correlation, and autocorrelation.

再次强调，卷积定理指出：

Again, the Convolution Theorem states:

计算卷积的一种方法是：

So one way to compute a convolution is:

其中 IDFT 是逆 DFT。卷积的简单实现所需时间与N²成正比；该算法使用 FFT，所需时间与 N²^成正比。

where IDFT is the inverse DFT. A simple implementation of convolution takes time proportional to N²; this algorithm, using the FFT, takes time proportional to .

我们可以通过双向计算相同的卷积来验证其有效性。例如，我将把它应用于图 8-1所示的 Facebook 股票数据：

We can confirm that it works by computing the same convolution both ways. As an example, I’ll apply it to the Facebook stock data shown in Figure 8-1:

导入 pandas as pd

names = ['日期', '开盘价', '最高价', '最低价', '收盘价', '成交量']
df = pd.read_csv('fb.csv', header=0, names=names)
ys = df.close.values[::-1]

import pandas as pd

names = ['date', 'open', 'high', 'low', 'close', 'volume']
df = pd.read_csv('fb.csv', header=0, names=names)
ys = df.close.values[::-1]

本示例使用 Pandas 从 CSV 文件（包含在本书的存储库中）读取数据。如果您不熟悉 Pandas，请不要担心：本书中我不会过多地使用它。但如果您感兴趣，可以在Think Stats （http://thinkstats2.com）上了解更多信息。

This example uses Pandas to read the data from the CSV file (included in the repository for this book). If you are not familiar with Pandas, don’t worry: I’m not going to do much with it in this book. But if you’re interested, you can learn more about it in Think Stats at http://thinkstats2.com.

结果df是一个 NumPy 数组，它是DataFramePandas 提供的数据结构之一，包含每日收盘价。close

The result, df, is a DataFrame, one of the data structures provided by Pandas. close is a NumPy array that contains daily closing prices.

接下来，我将创建一个高斯窗口，并将其与以下函数进行卷积close：

Next I’ll create a Gaussian window and convolve it with close:

window = scipy.signal.gaussian(M=30, std=6)
window /= window.sum()
smoothed = np.convolve(ys, window, mode='valid')

window = scipy.signal.gaussian(M=30, std=6)
window /= window.sum()
smoothed = np.convolve(ys, window, mode='valid')

fft_convolve使用快速傅里叶变换 (FFT) 计算相同的结果：

fft_convolve computes the same thing using the FFT:

from np.fft import fft, ifft

def fft_convolve(signal, window):
    fft_signal = fft(signal)
    fft_window = fft(window)
    返回 ifft(fft_signal * fft_window)

from np.fft import fft, ifft

def fft_convolve(signal, window):
    fft_signal = fft(signal)
    fft_window = fft(window)
    return ifft(fft_signal * fft_window)

我们可以通过将窗口填充到与原来的长度相同的长度ys，然后计算卷积来测试它：

We can test it by padding the window to the same length as ys and then computing the convolution:

padded = zero_pad(window, N)
smoothed2 = fft_convolve(ys, padded)

padded = zero_pad(window, N)
smoothed2 = fft_convolve(ys, padded)

结果开头部分存在错误值，其中M为窗口长度。我们可以像这样去除这些错误值：

The result has bogus values at the beginning, where M is the length of the window. We can slice off the bogus values like this:

M = len(window)
smoothed2 = smoothed2[M-1:]

M = len(window)
smoothed2 = smoothed2[M-1:]

结果与预期相符，fft_convolve精度约为 12 位有效数字。

The result agrees with fft_convolve with about 12 digits of precision.

高效自相关

Efficient Autocorrelation

在“卷积”一节中，我介绍了互相关和卷积的定义，我们看到它们几乎相同，只是卷积中的窗口是反转的。

In “Convolution” I presented definitions of cross-correlation and convolution, and we saw that they are almost the same, except that in convolution the window is reversed.

现在我们有了高效的卷积算法，也可以用它来计算互相关和自相关。利用上一节的数据，我们可以计算Facebook股票价格的自相关性：

Now that we have an efficient algorithm for convolution, we can also use it to compute cross-correlations and autocorrelations. Using the data from the previous section, we can compute the autocorrelation of Facebook stock prices:

corrs = np.correlate(close, close, mode='same')

corrs = np.correlate(close, close, mode='same')

当 时mode='same'，结果与 的长度相同close，对应于从到 的滞后。图 8-8中的灰线显示了该结果。除了 处之外，没有峰值，因此该信号没有明显的周期性行为。然而，自相关函数缓慢下降，表明该信号类似于粉红噪声，正如我们在“自相关”部分所看到的。lag=0

With mode='same', the result has the same length as close, corresponding to lags from to . The gray line in Figure 8-8 shows the result. Except at lag=0, there are no peaks, so there is no apparent periodic behavior in this signal. However, the autocorrelation function drops off slowly, suggesting that this signal resembles pink noise, as we saw in “Autocorrelation”.

图 8-8.由 NumPy 和 fft_autocorr 计算的自相关函数。

要使用卷积计算自相关，我们需要对信号进行零填充，使其长度翻倍。这是必要的，因为快速傅里叶变换 (FFT) 基于信号周期性的假设，即信号会从末尾循环到开头。对于像这样的时间序列数据，这个假设并不成立。通过添加零并截断结果，可以去除这些无效值。

To compute autocorrelation using convolution, we have to zero-pad the signal to double the length. This trick is necessary because the FFT is based on the assumption that the signal is periodic; that is, that it wraps around from the end to the beginning. With time-series data like this, that assumption is invalid. Adding zeros, and then trimming the results, removes the bogus values.

另外，请记住卷积会反转窗口方向。为了消除这种影响，我们在调用 `flash()` 之前先反转窗口方向fft_convolve，使用np.flipud`flash()` 函数翻转 NumPy 数组。结果是数组的视图，而不是副本，因此此操作速度很快：

Also, remember that convolution reverses the direction of the window. In order to cancel that effect, we reverse the direction of the window before calling fft_convolve, using np.flipud, which flips a NumPy array. The result is a view of the array, not a copy, so this operation is fast:

def fft_autocorr(signal):
    N = len(信号)
    signal = thinkdsp.zero_pad(signal, 2*N)
    window = np.flipud(signal)

    corrs = fft_convolve(signal, window)
    corrs = np.roll(corrs, N//2+1)[:N]
    返回值

def fft_autocorr(signal):
    N = len(signal)
    signal = thinkdsp.zero_pad(signal, 2*N)
    window = np.flipud(signal)

    corrs = fft_convolve(signal, window)
    corrs = np.roll(corrs, N//2+1)[:N]
    return corrs

结果fft_convolve长度为2N。其中，第一个和最后一个元素是有效的；其余元素是零填充的结果。为了选择有效元素，我们对结果进行滚动，并选择前N 个元素，对应于从到的滞后值。

The result from fft_convolve has length 2N. Of those, the first and last are valid; the rest are the result of zero-padding. To select the valid element, we roll the results and select the first N, corresponding to lags from to .

如图 8-8所示，fft_autocorr和 的结果np.correlate相同（精度约为 9 位有效数字）。

As shown in Figure 8-8, the results from fft_autocorr and np.correlate are identical (with about 9 digits of precision).

请注意，图 8-8中的相关性数值很大；我们可以按照“使用 NumPy”中所示，将它们归一化（介于 -1 和 1 之间） 。

Notice that the correlations in Figure 8-8 are large numbers; we could normalize them (between –1 and 1) as shown in “Using NumPy”.

我们这里用于自相关的策略也适用于互相关。同样，你需要先对信号进行预处理，方法是翻转其中一个信号并对两个信号都进行填充，然后再去除结果中的无效部分。这种填充和去除操作比较繁琐，但正因如此，像 NumPy 这样的库才提供了相应的函数来帮你完成这些工作。

The strategy we used here for autocorrelation also works for cross-correlation. Again, you have to prepare the signals by flipping one and padding both, and then you have to trim the invalid parts of the result. This padding and trimming is a nuisance, but that’s why libraries like NumPy provide functions to do it for you.

练习

Exercises

这些练习的答案在chap08soln.ipynb……

Solutions to these exercises are in chap08soln.ipynb.

有限差分法

Finite Differences

在“平滑”部分，我们对 Facebook 的每日股价应用了平滑窗口，发现时域中的平滑窗口对应于频域中的低通滤波器。

In “Smoothing”, we applied a smoothing window to the daily stock price of Facebook and found that a smoothing window in the time domain corresponds to a low-pass filter in the frequency domain.

在本节中，我们将观察每日价格变化，并发现计算时域中连续元素之间的差异相当于高通滤波器。

In this section, we’ll look at daily price changes and see that computing the difference between successive elements, in the time domain, corresponds to a high-pass filter.

以下是读取数据、将其存储为波形并计算其频谱的代码：

Here’s the code to read the data, store it as a wave, and compute its spectrum:

导入 pandas as pd

names = ['日期', '开盘价', '最高价', '最低价', '收盘价', '成交量']
df = pd.read_csv('fb.csv', header=0, names=names)
ys = df.close.values[::-1]
close = thinkdsp.Wave(ys, framerate=1)
spectrum = wave.make_spectrum()

import pandas as pd

names = ['date', 'open', 'high', 'low', 'close', 'volume']
df = pd.read_csv('fb.csv', header=0, names=names)
ys = df.close.values[::-1]
close = thinkdsp.Wave(ys, framerate=1)
spectrum = wave.make_spectrum()

本示例使用 Pandas 读取 CSV 文件；结果是一个DataFrame包含df开盘价、收盘价、最高价和最低价列的表格。我选择收盘价并将其保存到一个Wave对象中。帧速率为每天 1 个样本。

This example uses Pandas to read the CSV file; the result is a DataFrame, df, with columns for the opening price, closing price, and high and low prices. I select the closing prices and save them in a Wave object. The frame rate is 1 sample per day.

图 9-1显示了该时间序列及其频谱。从视觉上看，该时间序列类似于布朗噪声（参见“布朗噪声”）。频谱看起来像一条直线，尽管存在噪声。估计的斜率为 -1.9，这与布朗噪声的特征相符。

Figure 9-1 shows this time series and its spectrum. Visually, the time series resembles Brownian noise (see “Brownian Noise”). And the spectrum looks like a straight line, albeit a noisy one. The estimated slope is –1.9, which is consistent with Brownian noise.

图 9-1. Facebook 股票的每日收盘价及该时间序列的频谱。

现在让我们用以下公式计算每日价格变化np.diff：

Now let’s compute the daily price change using np.diff:

diff = np.diff(ys)
change = thinkdsp.Wave(diff, framerate=1)
change_spectrum = change.make_spectrum()

diff = np.diff(ys)
change = thinkdsp.Wave(diff, framerate=1)
change_spectrum = change.make_spectrum()

图 9-2显示了由此产生的波形及其频谱。每日变化类似于白噪声，频谱的估计斜率 -0.06 接近于零，这与我们对白噪声的预期相符。

Figure 9-2 shows the resulting wave and its spectrum. The daily changes resemble white noise, and the estimated slope of the spectrum, –0.06, is near zero, which is what we expect for white noise.

图 9-2. Facebook 股票每日价格变化及该时间序列的频谱。

频域

The Frequency Domain

计算相邻元素之间的差值与使用窗口进行卷积相同[1, -1]。如果这些元素的顺序看起来是反的，请记住，卷积运算在应用于信号之前会反转窗口。

Computing the difference between successive elements is the same as convolution with the window [1, -1]. If the order of those elements seems backward, remember that convolution reverses the window before applying it to the signal.

我们可以通过计算窗口的离散傅里叶变换（DFT）来观察该操作在频域中的效果：

We can see the effect of this operation in the frequency domain by computing the DFT of the window:

diff_window = np.array([1.0, -1.0])
padded = thinkdsp.zero_pad(diff_window, len(close))
diff_wave = thinkdsp.Wave(padded, framerate=close.framerate)
diff_filter = diff_wave.make_spectrum()

diff_window = np.array([1.0, -1.0])
padded = thinkdsp.zero_pad(diff_window, len(close))
diff_wave = thinkdsp.Wave(padded, framerate=close.framerate)
diff_filter = diff_wave.make_spectrum()

图 9-3显示了结果。有限差分窗口对应于一个高通滤波器：其幅度随频率增加而增加，低频时呈线性增加，高频后呈亚线性增加。下一节我们将解释原因。

Figure 9-3 shows the result. The finite difference window corresponds to a high-pass filter: its amplitude increases with frequency, linearly for low frequencies and then sublinearly after that. In the next section, we’ll see why.

图 9-3.对应于 diff 和 differential 算子（左）和积分算子（右，log-y 刻度）的滤波器。

分化

Differentiation

我们在上一节中使用的窗口是一阶导数的数值近似值，因此该滤波器近似于微分的效果。

The window we used in the previous section is a numerical approximation of the first derivative, so the filter approximates the effect of differentiation.

时域中的微分对应于频域中的简单滤波器；我们可以通过一些数学运算来弄清楚它是什么。

Differentiation in the time domain corresponds to a simple filter in the frequency domain; we can figure out what it is with a little math.

假设我们有一个频率为f 的复正弦波：

Suppose we have a complex sinusoid with frequency f:

E _f的一阶导数为：

The first derivative of E_f is:

我们可以将其改写为：

which we can rewrite as:

换句话说，对E _f求导就等于乘以，这是一个具有模和角度的复数。

In other words, taking the derivative of E_f is the same as multiplying by , which is a complex number with magnitude and angle .

我们可以这样计算与微分对应的滤波器：

We can compute the filter that corresponds to differentiation like this:

deriv_filter = close.make_spectrum()
deriv_filter.hs = PI2 * 1j * deriv_filter.fs

deriv_filter = close.make_spectrum()
deriv_filter.hs = PI2 * 1j * deriv_filter.fs

close我首先得到了具有合适尺寸和帧速率的频谱，然后hs用替换了。图 9-3（左）显示了该滤波器；它是一条直线。

I started with the spectrum of close, which has the right size and frame rate, then replaced the hs with . Figure 9-3 (left) shows this filter; it is a straight line.

正如我们在“矩阵合成”中看到的那样，将复正弦波乘以复数会产生两个效果：它会将振幅乘以 ，在本例中乘以，并将相位偏移移动，在本例中移动。

As we saw in “Synthesis with Matrices”, multiplying a complex sinusoid by a complex number has two effects: it multiplies the amplitude, in this case by , and shifts the phase offset, in this case by .

如果您熟悉算子和特征函数的语言，那么每个E _f都是微分算子的特征函数，其对应的特征值为。参见http://en.wikipedia.org/wiki/Eigenfunction。

If you are familiar with the language of operators and eigenfunctions, each E_f is an eigenfunction of the differentiation operator, with the corresponding eigenvalue . See http://en.wikipedia.org/wiki/Eigenfunction.

如果您不熟悉这种语言，它的意思如下：

If you are not familiar with that language, here’s what it means:

• 运算符是一种接受一个函数并返回另一个函数的函数。例如，微分就是一个运算符。

• An operator is a function that takes a function and returns another function. For example, differentiation is an operator.

• 如果对函数g应用算子 会产生函数乘以一个标量的效果，则称函数g是算子 的本征函数。也就是说，。

• A function, g, is an eigenfunction of an operator, , if applying to g has the effect of multiplying the function by a scalar. That is, .

• 在这种情况下，标量λ是对应于特征函数g 的特征值。

• In that case, the scalar λ is the eigenvalue that corresponds to the eigenfunction g.

• 一个给定的算符可能有很多本征函数，每个本征函数都有一个对应的本征值。

• A given operator might have many eigenfunctions, each with a corresponding eigenvalue.

因为复正弦函数是微分算子的特征函数，所以很容易求导。我们只需要乘以一个复标量即可。

Because complex sinusoids are eigenfunctions of the differentiation operator, they are easy to differentiate. All we have to do is multiply by a complex scalar.

对于包含多个分量的信号，处理过程只是稍微复杂一些：

For signals with more than one component, the process is only slightly harder:

1. 将信号表示为复正弦波之和。

2. Express the signal as the sum of complex sinusoids.

3. 通过乘法计算每个分量的导数。

4. Compute the derivative of each component by multiplication.

5. 将各差异化组成部分相加。

6. Add up the differentiated components.

如果这个过程听起来很熟悉，那是因为它与“高效卷积”中的卷积算法完全相同：计算 DFT，乘以滤波器，然后计算逆 DFT。

If that process sounds familiar, that’s because it is identical to the algorithm for convolution in “Efficient Convolution”: compute the DFT, multiply by a filter, and compute the inverse DFT.

Spectrum提供了一种应用微分滤波器的方法：

Spectrum provides a method that applies the differentiation filter:

# 光谱类：

    def different(self):
        self.hs *= PI2 * 1j * self.fs

# class Spectrum:

    def differentiate(self):
        self.hs *= PI2 * 1j * self.fs

我们可以用它来计算Facebook时间序列的导数：

We can use it to compute the derivative of the Facebook time series:

deriv_spectrum = close.make_spectrum()
deriv_spectrum.differentiate()
deriv = deriv_spectrum.make_wave()

deriv_spectrum = close.make_spectrum()
deriv_spectrum.differentiate()
deriv = deriv_spectrum.make_wave()

图 9-4将通过计算得到的每日价格变化np.diff与我们刚刚计算的导数进行了比较。我选择了时间序列中的前 50 个值，以便更清楚地看到差异。

Figure 9-4 compares the daily price changes computed by np.diff with the derivative we just computed. I selected the first 50 values in the time series so we can see the differences more clearly.

图 9-4. np.diff 计算的每日价格变化与应用微分滤波器计算的每日价格变化的比较。

如图 9-3 （左）所示，导数噪声较大，因为它放大了高频分量。此外，导数的前几个元素噪声也很大。问题在于，基于 DFT 的导数是基于信号周期性的假设。实际上，它将时间序列中的最后一个元素与第一个元素连接起来，从而在边界处产生伪影。

The derivative is noisier, because it amplifies the high-frequency components more, as shown in Figure 9-3 (left). Also, the first few elements of the derivative are very noisy. The problem there is that the DFT-based derivative is based on the assumption that the signal is periodic. In effect, it connects the last element in the time series back to the first element, which creates artifacts at the boundaries.

综上所述，我们已经证明：

To summarize, we have shown:

• 计算信号中连续值之间的差值可以表示为与简单窗口进行卷积。结果是对一阶导数的近似值。

• Computing the difference between successive values in a signal can be expressed as convolution with a simple window. The result is an approximation of the first derivative.

• 时域微分对应于频域中的简单滤波器。对于周期信号，其结果恰好是其一阶导数。对于某些非周期信号，它可以近似地表示导数。

• Differentiation in the time domain corresponds to a simple filter in the frequency domain. For periodic signals, the result is the first derivative, exactly. For some non-periodic signals, it can approximate the derivative.

利用 DFT 计算导数是求解微分方程的谱方法的基础（参见http://en.wikipedia.org/wiki/Spectral_method）。

Using the DFT to compute derivatives is the basis of spectral methods for solving differential equations (see http://en.wikipedia.org/wiki/Spectral_method).

特别是，它对分析线性时不变系统很有用，这将在第 10 章中讨论。

In particular, it is useful for the analysis of linear, time-invariant systems, which is coming up in Chapter 10.

一体化

Integration

在前一节中，我们证明了时域中的微分对应于频域中的一个简单的滤波器：它将每个分量乘以。由于积分是微分的逆运算，它也对应于一个简单的滤波器：它将每个分量除以。

In the previous section, we showed that differentiation in the time domain corresponds to a simple filter in the frequency domain: it multiplies each component by . Since integration is the inverse of differentiation, it also corresponds to a simple filter: it divides each component by .

我们可以这样计算这个滤波器：

We can compute this filter like this:

integ_filter = close.make_spectrum()
integ_filter.hs = 1 / (PI2 * 1j * integ_filter.fs)

integ_filter = close.make_spectrum()
integ_filter.hs = 1 / (PI2 * 1j * integ_filter.fs)

图 9-3（右）以对数y刻度显示了该滤波器，这样更容易看清。

Figure 9-3 (right) shows this filter on a log-y scale, which makes it easier to see.

Spectrum提供了一种应用积分过滤器的方法：

Spectrum provides a method that applies the integration filter:

# 光谱类：

    def integrate(self):
        self.hs /= PI2 * 1j * self.fs

# class Spectrum:

    def integrate(self):
        self.hs /= PI2 * 1j * self.fs

我们可以通过将积分滤波器应用于我们刚刚计算出的导数的频谱来验证其正确性：

We can confirm that the integration filter is correct by applying it to the spectrum of the derivative we just computed:

integ_spectrum = deriv_spectrum.copy()
integ_spectrum.integrate()

integ_spectrum = deriv_spectrum.copy()
integ_spectrum.integrate()

但请注意，在除法运算中，我们除以了 0。在 NumPy 中，结果是 NaN，这是一个特殊的浮点值，表示“非数字”。我们可以通过在将频谱转换回波形之前将此值设置为 0 来部分解决这个问题：

But notice that at , we are dividing by 0. The result in NumPy is NaN, which is a special floating-point value that represents “not a number”. We can partially deal with this problem by setting this value to 0 before converting the spectrum back to a wave:

integ_spectrum.hs[0] = 0
integ_wave = integ_spectrum.make_wave()

integ_spectrum.hs[0] = 0
integ_wave = integ_spectrum.make_wave()

图 9-5显示了积分导数及其原始时间序列。它们几乎完全相同，但积分导数向下平移了。问题在于，当我们处理该分量时，我们将信号的偏置设置为 0。但这并不奇怪；通常，微分会丢失偏置信息，而积分无法恢复它。从某种意义上说，此处的 NaN 值表明该元素未知。

Figure 9-5 shows this integrated derivative along with the original time series. They are almost identical, but the integrated derivative has been shifted down. The problem is that when we clobbered the component, we set the bias of the signal to 0. But that should not be surprising; in general, differentiation loses information about the bias, and integration can’t recover it. In some sense, the NaN at is telling us that this element is unknown.

图 9-5.原始时间序列与积分导数的比较。

如果我们提供这个“积分常数”，结果就相同了，这证实了这个积分滤波器是微分滤波器的正确逆滤波器。

If we provide this “constant of integration”, the results are identical, which confirms that this integration filter is the correct inverse of the differentiation filter.

累计总和

Cumulative Sum

就像微分算子近似于积分一样diff，累积和近似于积分。我将用锯齿波信号来演示：

In the same way that the diff operator approximates differentiation, the cumulative sum approximates integration. I’ll demonstrate with a sawtooth signal:

signal = thinkdsp.SawtoothSignal(freq=50)
in_wave = signal.make_wave(duration=0.1, framerate=44100)

signal = thinkdsp.SawtoothSignal(freq=50)
in_wave = signal.make_wave(duration=0.1, framerate=44100)

图 9-6显示了该波及其频谱。

Figure 9-6 shows this wave and its spectrum.

Wave提供一个计算波形数组累积和并返回新Wave对象的方法：

Wave provides a method that computes the cumulative sum of a wave array and returns a new Wave object:

# Wave 类：

    def cumsum(self):
        ys = np.cumsum(self.ys)
        ts = self.ts.copy()
        返回 Wave(ys, ts, self.framerate)

# class Wave:

    def cumsum(self):
        ys = np.cumsum(self.ys)
        ts = self.ts.copy()
        return Wave(ys, ts, self.framerate)

我们可以用它来计算以下各项的累积和in_wave：

We can use it to compute the cumulative sum of in_wave:

输出波 = 输入波.累积值()
out_wave.unbias()

out_wave = in_wave.cumsum()
out_wave.unbias()

图 9-7显示了所得波形及其频谱。如果您完成了第二章的练习，那么这个波形应该看起来很眼熟：它是一个抛物线信号。

Figure 9-7 shows the resulting wave and its spectrum. If you did the exercises in Chapter 2, this waveform should look familiar: it’s a parabolic signal.

图 9-6.锯齿波及其频谱。

图 9-7.抛物线波及其频谱。

比较两者，我们发现抛物线信号频谱中各分量的振幅衰减速度比锯齿波信号频谱中更快。在第二章中，我们看到锯齿波分量的衰减与 成正比。由于累积和近似于积分，而积分会按 成比例地滤波分量，因此抛物线波分量的衰减与 成正比。

Comparing the two, we see that the amplitudes of the components drop off more quickly in the spectrum of the parabolic signal than in the spectrum of the sawtooth signal. In Chapter 2, we saw that the components of the sawtooth drop off in proportion to . Since the cumulative sum approximates integration, and integration filters components in proportion to , the components of the parabolic wave drop off in proportion to .

我们可以通过计算与累积和对应的滤波器，以图形方式看到这一点：

We can see that graphically by computing the filter that corresponds to the cumulative sum:

cumsum_filter = diff_filter.copy()
cumsum_filter.hs = 1 / cumsum_filter.hs

cumsum_filter = diff_filter.copy()
cumsum_filter.hs = 1 / cumsum_filter.hs

因为cumsum是 的逆运算diff，所以我们从 的副本开始diff_filter，它是与运算 对应的滤波器diff，然后对 进行逆运算hs。

Because cumsum is the inverse operation of diff, we start with a copy of diff_filter, which is the filter that corresponds to the diff operation, and then invert the hs.

图 9-8显示了对应于累积和滤波器和积分滤波器的滤波器。累积和滤波器是对积分滤波器的良好近似，但在最高频率处，其衰减速度略快。

Figure 9-8 shows the filters corresponding to cumulative sum and integration. The cumulative sum is a good approximation of integration except at the highest frequencies, where it drops off a little faster.

图 9-8.与累积和及积分对应的滤波器。

为了确认这是累积和的正确滤波器，我们可以将其out_wave与以下频谱的比值进行比较in_wave：

To confirm that this is the correct filter for the cumulative sum, we can compare it to the ratio of the spectrum out_wave to the spectrum of in_wave:

in_spectrum = in_wave.make_spectrum()
out_spectrum = out_wave.make_spectrum()
ratio_spectrum = out_spectrum.ratio(in_spectrum, thresh=1)

in_spectrum = in_wave.make_spectrum()
out_spectrum = out_wave.make_spectrum()
ratio_spectrum = out_spectrum.ratio(in_spectrum, thresh=1)

以下是计算这些比率的方法：

And here’s the method that computes the ratios:

def ratio(self, denom, thresh=1):
    ratio_spectrum = self.copy()
    ratio_spectrum.hs /= denom.hs
    ratio_spectrum.hs[denom.amps < thresh] = np.nan
    返回比率光谱

def ratio(self, denom, thresh=1):
    ratio_spectrum = self.copy()
    ratio_spectrum.hs /= denom.hs
    ratio_spectrum.hs[denom.amps < thresh] = np.nan
    return ratio_spectrum

当denom.amps较小时，所得比率会有噪声，所以我将这些值设置为 NaN。

When denom.amps is small, the resulting ratio is noisy, so I set those values to NaN.

图 9-9显示了累积和对应的比率和滤波器。它们吻合，这证实了反转滤波器即可diff得到滤波器cumsum。

Figure 9-9 shows the ratios and the filter corresponding to the cumulative sum. They agree, which confirms that inverting the filter for diff yields the filter for cumsum.

图 9-9.与前后光谱的累积和及实际比率相对应的滤波器。

cumsum最后，我们可以通过在频域中应用滤波器来验证卷积定理的适用性：

Finally, we can confirm that the Convolution Theorem applies by applying the cumsum filter in the frequency domain:

out_wave2 = (in_spectrum * cumsum_filter).make_wave()

out_wave2 = (in_spectrum * cumsum_filter).make_wave()

在浮点误差范围内，它与我们用计算得到的out_wave2完全相同，因此卷积定理成立！但请注意，此演示仅适用于周期信号。out_wavecumsum

Within the limits of floating-point error, out_wave2 is identical to out_wave, which we computed using cumsum, so the Convolution Theorem works! But note that this demonstration only works with periodic signals.

噪声积分

Integrating Noise

在“布朗噪声”一节中，我们通过计算白噪声的累积和生成了布朗噪声。现在我们已经了解了cumsum频域中的影响，对布朗噪声的频谱也有了一定的认识。

In “Brownian Noise”, we generated Brownian noise by computing the cumulative sum of white noise. Now that we understand the effect of cumsum in the frequency domain, we have some insight into the spectrum of Brownian noise.

白噪声在所有频率上的平均功率相等。计算累积和时，每个分量的幅度除以频率f。由于功率是幅度的平方，因此每个分量的功率除以f² 。所以，平均而言，频率f处的功率与^以下公式成正比：

White noise has equal power at all frequencies, on average. When we compute the cumulative sum, the amplitude of each component is divided by f. Since power is the square of magnitude, the power of each component is divided by f². So on average, the power at frequency f is proportional to :

其中K是一个无关紧要的常数。两边取对数得到：

where K is a constant that’s not important. Taking the log of both sides yields:

因此，当我们在双对数坐标系中绘制布朗噪声的频谱时，我们期望看到一条斜率为 -2 的直线，至少近似如此。

And that’s why, when we plot the spectrum of Brownian noise on a log-log scale, we expect to see a straight line with slope –2, at least approximately.

在“有限差分法”中，我们绘制了Facebook股票收盘价的频谱图，并估计其斜率为-1.9，这与布朗噪声一致。许多股票价格都具有类似的频谱。

In “Finite Differences” we plotted the spectrum of closing prices for Facebook stock, and estimated that the slope is –1.9, which is consistent with Brownian noise. Many stock prices have similar spectrums.

当我们使用该diff算子计算每日变化时，我们将每个分量的振幅乘以一个与f成正比的滤波器，这意味着我们将每个分量的功率乘以f²。在双对数坐标系中，此操作将功率谱的斜率加 2，这就是为什么结果的估计斜率接近 0.1（但略低一些，因为它^只是diff近似微分）。

When we used the diff operator to compute daily changes, we multiplied the amplitude of each component by a filter proportional to f, which means we multiplied the power of each component by f². On a log-log scale, this operation adds 2 to the slope of the power spectrum, which is why the estimated slope of the result is near 0.1 (but a little lower, because diff only approximates differentiation).

练习

Exercises

这些练习的答案在chap09soln.ipynb……

Solutions to these exercises are in chap09soln.ipynb.

信号与系统

Signals and Systems

在信号处理的语境中，系统是对任何以信号作为输入并产生信号作为输出的事物的抽象表示。

In the context of signal processing, a system is an abstract representation of anything that takes a signal as input and produces a signal as output.

例如，电子放大器是一种以电信号作为输入并产生（更大的）信号作为输出的电路。

For example, an electronic amplifier is a circuit that takes an electrical signal as input and produces a (louder) signal as output.

再举一个例子，当你听音乐表演时，你可以把房间想象成一个系统，它会接收表演产生位置的声音，并在你听到声音的位置产生略有不同的声音。

As another example, when you listen to a musical performance, you can think of the room as a system that takes the sound of the performance at the location where it is generated and produces a somewhat different sound at the location where you hear it.

线性时不变系统¹是具有以下两个性质的系统：

A linear, time-invariant system¹ is a system with these two properties:

线性

Linearity

如果同时向系统中输入两个信号，则结果等于它们的输出之和。数学上，如果输入x₁产生输出_{y₁ ，}另一个输入x₂产生输出y₂，则x₁ + _y₂ = y₁ + _y₂，其中a和b是标量。

If you put two inputs into the system at the same time, the result is the sum of their outputs. Mathematically, if an input x₁ produces output y₁ and another input x₂ produces y₂, then produces , where a and b are scalars.

时间不变性

Time invariance

系统的效果不会随时间变化，也不取决于系统的状态。因此，如果输入_x1和 x2 仅在时间上有所不同，则它们的输出 y1 和 y2也相差相同的_时间差，_{但除此之外}完全相同。

The effect of the system doesn’t vary over time, or depend on the state of the system. So if inputs x₁ and x₂ differ only by a shift in time, their outputs, y₁ and y₂, differ by the same shift but are otherwise identical.

许多物理系统都具有这些特性，至少近似如此：

Many physical systems have these properties, at least approximately:

• 仅包含电阻器、电容器和电感器的电路是线性时不变 (LTI) 电路，因为这些元件的行为与其理想模型一致。

• Circuits that contain only resistors, capacitors, and inductors are LTI, to the degree that the components behave like their idealized models.

• 包含弹簧、质量块和阻尼器的机械系统也是线性时不变系统，假设弹簧是线性的（力与位移成正比），阻尼器是线性的（力与速度成正比）。

• Mechanical systems that contain springs, masses, and dashpots are also LTI, assuming linear springs (force proportional to displacement) and dashpots (force proportional to velocity).

• 此外，与本书中的应用最为相关的是，传播声音的介质（包括空气、水和固体）可以用 LTI 系统很好地建模。

• Also, and most relevant to applications in this book, the media that transmit sound (including air, water, and solids) are well modeled by LTI systems.

LTI 系统由线性微分方程描述，这些方程的解是复正弦波（参见http://en.wikipedia.org/wiki/Linear_differential_equation）。

LTI systems are described by linear differential equations, and the solutions of those equations are complex sinusoids (see http://en.wikipedia.org/wiki/Linear_differential_equation).

该结果提供了一种计算线性时不变系统对输入信号影响的算法：

This result provides an algorithm for computing the effect of an LTI system on an input signal:

1. 将信号表示为复正弦分量之和。

2. Express the signal as the sum of complex sinusoid components.

3. 对于每个输入分量，计算相应的输出分量。

4. For each input component, compute the corresponding output component.

5. 将输出分量相加。

6. Add up the output components.

到目前为止，我希望这个算法听起来很熟悉。它与我们在“高效卷积”一节中用于卷积的算法以及在“微分”一节中用于微分的算法相同。这个过程被称为谱分解，因为我们将输入信号“分解”成它的频谱分量。

At this point, I hope this algorithm sounds familiar. It’s the same algorithm we used for convolution in “Efficient Convolution”, and for differentiation in “Differentiation”. This process is called spectral decomposition because we “decompose” the input signal into its spectral components.

为了将此方法应用于线性时不变（LTI）系统，我们必须通过找出系统对输入信号各分量的影响来表征该系统。对于机械系统，有一种简单有效的方法：踢一下它，然后记录输出。

In order to apply this process to an LTI system, we have to characterize the system by finding its effect on each component of the input signal. For mechanical systems, it turns out that there is a simple and efficient way to do that: you kick it and record the output.

从技术上讲，这种“冲击”被称为脉冲，而输出被称为脉冲响​​应。你可能会好奇，单个脉冲如何能完整地表征一个系统。你可以通过计算脉冲的离散傅里叶变换 (DFT) 来找到答案。以下是一个脉冲位于 处的波阵列：

Technically, the “kick” is called an impulse and the output is called the impulse response. You might wonder how a single impulse can completely characterize a system. You can see the answer by computing the DFT of an impulse. Here’s a wave array with an impulse at :

impulse = np.zeros(8)
impulse[0] = 1
impulse_spectrum = np.fft.fft(impulse)

impulse = np.zeros(8)
impulse[0] = 1
impulse_spectrum = np.fft.fft(impulse)

这是波形阵列：

Here’s the wave array:

[ 1. 0. 0. 0. 0. 0. 0. 0.]

[ 1.  0.  0.  0.  0.  0.  0.  0.]

以下是它的光谱：

And here’s its spectrum:

[ 1.+0.j 1.+0.j 1.+0.j 1.+0.j 1.+0.j 1.+0.j 1.+0.j 1.+0.j]

[ 1.+0.j  1.+0.j  1.+0.j  1.+0.j  1.+0.j  1.+0.j  1.+0.j  1.+0.j]

该频谱全为1；也就是说，一个脉冲是由所有频率上幅度相等的分量叠加而成的。这种频谱不应与白噪声混淆，白噪声在所有频率上的平均功率相同，但其值围绕该平均值波动。

The spectrum is all ones; that is, an impulse is the sum of components with equal magnitudes at all frequencies. This spectrum should not be confused with white noise, which has the same average power at all frequencies, but varies around that average.

当你通过输入脉冲信号来测试系统时，你实际上是在测试系统在所有频率下的响应。而且你可以同时测试所有频率，因为系统是线性的，所以同时进行的测试不会相互干扰。

When you test a system by inputting an impulse, you are testing the response of the system at all frequencies. And you can test them all at the same time because the system is linear, so simultaneous tests don’t interfere with each other.

窗口和过滤器

Windows and Filters

为了说明这种系统表征方法的有效性，我将从一个简单的例子开始：双元素移动平均。我们可以将这种运算视为一个系统，它以一个信号作为输入，并产生一个略微平滑的信号作为输出。

To show why this kind of system characterization works, I will start with a simple example: a 2-element moving average. We can think of this operation as a system that takes a signal as an input and produces a slightly smoother signal as an output.

在这个例子中，我们知道窗口大小，所以可以计算出相应的滤波器。但通常情况并非如此；下一节我们将看一个事先不知道窗口大小和滤波器大小的例子。

In this example we know what the window is, so we can compute the corresponding filter. But that’s not usually the case; in the next section we’ll look at an example where we don’t know the window or the filter ahead of time.

这里有一个计算双元素移动平均值的窗口（参见“平滑”）：

Here’s a window that computes a 2-element moving average (see “Smoothing”):

window_array = np.array([0.5, 0.5, 0, 0, 0, 0, 0, 0,])
window = thinkdsp.Wave(window_array, framerate=8)

window_array = np.array([0.5, 0.5, 0, 0, 0, 0, 0, 0,])
window = thinkdsp.Wave(window_array, framerate=8)

我们可以通过计算窗口的DFT来找到相应的滤波器：

We can find the corresponding filter by computing the DFT of the window:

filtr = window.make_spectrum(full=True)

filtr = window.make_spectrum(full=True)

图 10-1显示了结果。与移动平均窗口对应的滤波器是一个近似高斯曲线形状的低通滤波器。

Figure 10-1 shows the result. The filter that corresponds to a moving average window is a low-pass filter with the approximate shape of a Gaussian curve.

图 10-1. 2 元素移动平均窗口的 DFT。

现在假设我们不知道窗口函数或相应的滤波器，而我们想要表征这个系统。我们可以通过输入一个脉冲并测量脉冲响应来实现这一点。

Now imagine that we did not know the window or the corresponding filter, and we wanted to characterize this system. We would do that by inputting an impulse and measuring the impulse response.

在这个例子中，我们可以通过将脉冲的频谱与滤波器相乘来计算脉冲响应，然后将结果从频谱转换为波形：

In this example, we can compute the impulse response by multiplying the spectrum of the impulse and the filter, and then converting the result from a spectrum to a wave:

产物 = 脉冲频谱 * 滤波器
过滤后的 = product.make_wave()

product = impulse_spectrum * filtr
filtered = product.make_wave()

由于impulse_spectrum所有元素均为 1，因此该乘积与滤波器相同，而滤波后的波形与窗口相同。

Since impulse_spectrum is all ones, the product is identical to the filter, and the filtered wave is identical to the window.

这个例子说明了两件事：

This example demonstrates two things:

• 因为脉冲的频谱全为 1，所以脉冲响应的 DFT 与表征该系统的滤波器相同。

• Because the spectrum of an impulse is all ones, the DFT of the impulse response is identical to the filter that characterizes the system.

• 因此，脉冲响应与表征该系统的卷积窗口相同。

• Therefore, the impulse response is identical to the convolution window that characterizes the system.

声学响应

Acoustic Response

为了描述房间或开放空间的声学响应，一种简单的产生脉冲的方法是戳破气球或开枪。这样得到的输入信号近似于脉冲，因此你听到的声音也近似于脉冲响应。

To characterize the acoustic response of a room or open space, a simple way to generate an impulse is to pop a balloon or fire a gun. The result is an input signal that approximates an impulse, so the sound you hear approximates the impulse response.

例如，我将使用一段枪声录音来描述开枪房间的特征，然后使用脉冲响应来模拟该房间对小提琴录音的影响。

As an example, I’ll use a recording of a gunshot to characterize the room where the gun was fired, then use the impulse response to simulate the effect of that room on a violin recording.

此示例位于本书的存储库中；您也可以在http://tinyurl.com/thinkdsp10chap10.ipynb查看并收听示例。

This example is in chap10.ipynb, which is in the repository for this book; you can also view it, and listen to the examples, at http://tinyurl.com/thinkdsp10.

这是枪声：

Here’s the gunshot:

response = thinkdsp.read_wave('180961__kleeb__gunshots.wav')
response = response.segment(start=0.26, duration=5.0)
response.normalize()
response.plot()

response = thinkdsp.read_wave('180961__kleeb__gunshots.wav')
response = response.segment(start=0.26, duration=5.0)
response.normalize()
response.plot()

我选择从 0.26 秒开始的片段，去除枪声之前的静音部分。图 10-2（左）显示了枪声的波形。接下来，我计算以下波形的离散傅里叶变换 (DFT) response：

I select a segment starting at 0.26 seconds to remove the silence before the gunshot. Figure 10-2 (left) shows the waveform of the gunshot. Next I compute the DFT of response:

传输 = response.make_spectrum()
transfer.plot()

transfer = response.make_spectrum()
transfer.plot()

图 10-2.枪声波形。

图 10-2（右）显示了结果。该频谱编码了房间的响应；对于每个频率，频谱包含一个复数，分别表示幅度倍增器和相移。该频谱被称为传递函数，因为它包含了系统如何将输入传递到输出的信息。

Figure 10-2 (right) shows the result. This spectrum encodes the response of the room; for each frequency, the spectrum contains a complex number that represents an amplitude multiplier and a phase shift. This spectrum is called a transfer function because it contains information about how the system transfers the input to the output.

现在我们可以模拟这个房间对小提琴音色的影响。这是我们在“周期信号”课程中使用的小提琴录音。：

Now we can simulate the effect this room would have on the sound of a violin. Here is the violin recording we used in “Periodic Signals”:

violin = thinkdsp.read_wave('92002__jcveliz__violin-origional.wav')
violin.truncate(len(response))
小提琴.normalize()

violin = thinkdsp.read_wave('92002__jcveliz__violin-origional.wav')
violin.truncate(len(response))
violin.normalize()

小提琴声和枪声的采样帧率相同，均为 44,100 Hz。巧合的是，两者的持续时间也大致相同。我将小提琴声裁剪到与枪声相同的长度。

The violin and gunshot waves were sampled at the same frame rate, 44,100 Hz. And coincidentally, the duration of both is about the same. I trimmed the violin wave to the same length as the gunshot.

接下来，我计算小提琴波形的离散傅里叶变换：

Next I compute the DFT of the violin wave:

spectrum = violin.make_spectrum()

spectrum = violin.make_spectrum()

现在我知道了输入信号中每个频率分量的幅值和相位，也知道了系统的传递函数。它们的乘积就是输出信号的离散傅里叶变换（DFT），我们可以用它来计算输出波形：

Now I know the magnitude and phase of each frequency component in the input, and I know the transfer function of the system. Their product is the DFT of the output, which we can use to compute the output wave:

输出 = (频谱 * 传输).make_wave()
output.normalize()
output.plot()

output = (spectrum * transfer).make_wave()
output.normalize()
output.plot()

图 10-3显示了系统的输入（上）和输出（下）。它们之间存在显著差异，而且这些差异在听觉上非常明显。请加载chap10.ipynb并聆听。我发现这个例子最引人注目的一点是，你可以感受到房间的声学特性；在我听来，它就像一个狭长的房间，地板和天花板都很硬。也就是说，就像一个射击场。

Figure 10-3 shows the input (top) and output (bottom) of the system. They are substantially different, and the differences are clearly audible. Load chap10.ipynb and listen to them. One thing I find striking about this example is that you can get a sense of what the room is like; to me, it sounds like a long, narrow room with hard floors and ceilings. That is, like a firing range.

图 10-3.卷积前后小提琴录音的波形。

在这个例子中，我忽略了一个细节，为了避免有人注意到，我在这里提一下。我最初使用的小提琴录音已经被一个系统（即录音房间）转换过了。所以，我实际计算的是经过两次转换后的小提琴声音。为了更准确地模拟小提琴在不同房间里的声音，我应该先对录音房间进行特性分析，然后应用该传递函数的逆函数。

There’s one thing I glossed over in this example that I’ll mention in case it bothers anyone. The violin recording I started with has already been transformed by one system: the room where it was recorded. So what I really computed in my example is the sound of the violin after two transformations. To properly simulate the sound of a violin in a different room, I should have characterized the room where the violin was recorded and applied the inverse of that transfer function first.

系统与卷积

Systems and Convolution

如果你觉得前面的例子是黑魔法，你并不孤单。我思考这个问题很久了，至今仍然觉得头疼。

If you think the previous example is black magic, you are not alone. I’ve been thinking about it for a while and it still makes my head hurt.

在前一节中，我提出了一种思考方式：

In the previous section, I suggested one way to think about it:

• 脉冲由所有频率上振幅均为 1 的分量组成。

• An impulse is made up of components with amplitude 1 at all frequencies.

• 脉冲响应包含系统对所有这些分量的响应之和。

• The impulse response contains the sum of the responses of the system to all of these components.

• 传递函数（即脉冲响应的 DFT）以幅度倍增器和相移的形式编码了系统对每个频率分量的影响。

• The transfer function, which is the DFT of the impulse response, encodes the effect of the system on each frequency component in the form of an amplitude multiplier and a phase shift.

• 对于任何输入，我们可以通过将输入分解成各个组成部分，计算每个组成部分的响应，然后将它们相加来计算系统的响应。

• For any input, we can compute the response of the system by breaking the input into components, computing the response to each component, and adding them up.

但如果你不喜欢这种方式，还有另一种完全不同的思考方式：卷积！根据卷积定理，频域中的乘法对应于时域中的卷积。在这个例子中，系统的输出是输入和系统响应的卷积。

But if you don’t like that, there’s another way to think about it altogether: convolution! By the Convolution Theorem, multiplication in the frequency domain corresponds to convolution in the time domain. In this example, the output of the system is the convolution of the input and the system response.

以下是理解其原理的关键：

Here are the keys to understanding why that works:

• 你可以把输入波形中的样本看作是一系列振幅变化的脉冲。

• You can think of the samples in the input wave as a sequence of impulses with varying amplitude.

• 输入中的每个脉冲都会产生一个脉冲响应的副本，该副本在时间上发生偏移（因为系统是时不变的），并按输入的幅度进行缩放。

• Each impulse in the input yields a copy of the impulse response, shifted in time (because the system is time-invariant) and scaled by the amplitude of the input.

• 输出是脉冲响应的平移和缩放副本之和。这些副本相加是因为系统是线性的。

• The output is the sum of the shifted, scaled copies of the impulse response. The copies add up because the system is linear.

让我们循序渐进地进行。假设我们不发射一门枪，而是发射两门：一门是振幅为 1 的大枪，另一门是振幅为 0.5 的小枪。

Let’s work our way up gradually. Suppose that instead of firing one gun, we fire two: a big one with amplitude 1 at and a smaller one with amplitude 0.5 at .

我们可以通过将原始脉冲响应与其缩放和平移后的副本相加来计算系统的响应。以下是一个生成波形平移缩放副本的函数：

We can compute the response of the system by adding up the original impulse response and a scaled, shifted copy of itself. Here’s a function that makes a shifted, scaled copy of a wave:

def shifted_scaled(wave, shift, factor):
    res = wave.copy()
    res.shift(shift)
    res.scale(因子)
    返回结果

def shifted_scaled(wave, shift, factor):
    res = wave.copy()
    res.shift(shift)
    res.scale(factor)
    return res

该参数shift为时间偏移量（以秒为单位）；factor是一个乘法因子。

The parameter shift is a time shift in seconds; factor is a multiplicative factor.

以下是我们如何使用它来计算对两响礼炮的回应：

Here’s how we use it to compute the response to a two-gun salute:

偏移量 = 1
系数 = 0.5
gun2 = response + shifted_scaled(response, shift, factor)

shift = 1
factor = 0.5
gun2 = response + shifted_scaled(response, shift, factor)

图 10-4显示了结果。你可以听听它的声音chap10.ipynb。不出所料，听起来像是两声枪响，第一声比第二声响。

Figure 10-4 shows the result. You can hear what it sounds like in chap10.ipynb. Not surprisingly, it sounds like two gunshots, the first one louder than the second.

图 10-4.波形与平移、缩放副本之和。

现在假设你不是用两支枪，而是用一百支枪，以每秒 441 发的速率进行射击。这个循环计算结果：

Now suppose instead of 2 guns, you were to add up 100 guns fired at a rate of 441 shots per second. This loop computes the result:

dt = 1 / 441
总计 = 0
for k in range(100):
    total += shifted_scaled(response, k*dt, 1.0)

dt = 1 / 441
total = 0
for k in range(100):
    total += shifted_scaled(response, k*dt, 1.0)

每秒发射 441 次，你听不到每次射击的声音。取而代之的是，它听起来像是一个频率为 441 赫兹的周期性信号。如果你播放这个例子，它听起来就像车库里的汽车喇叭声。

With 441 shots per second, you don’t hear the individual shots. Instead, it sounds like a periodic signal at 441 Hz. If you play this example, it sounds like a car horn in a garage.

这就引出了一个关键见解：你可以把任何波看作一系列样本，其中每个样本都是一个振幅不同的脉冲。

And that brings us to a key insight: you can think of any wave as a series of samples, where each sample is an impulse with a different amplitude.

例如，我将生成一个频率为 441 Hz 的锯齿波信号：

As an example, I’ll generate a sawtooth signal at 441 Hz:

signal = thinkdsp.SawtoothSignal(freq=441)
wave = signal.make_wave(duration=0.1,
                        帧率=响应.帧率)

signal = thinkdsp.SawtoothSignal(freq=441)
wave = signal.make_wave(duration=0.1,
                        framerate=response.framerate)

现在我将遍历构成锯齿波的一系列脉冲，并将脉冲响应相加：

Now I’ll loop through the series of impulses that make up the sawtooth, and add up the impulse responses:

总计 = 0
for t, y in zip(wave.ts, wave.ys):
    total += shifted_scaled(response, t, y)

total = 0
for t, y in zip(wave.ts, wave.ys):
    total += shifted_scaled(response, t, y)

结果听起来就像在射击场上播放锯齿波一样。你可以再次收听chap10.ipynb。

The result is what it would sound like to play a sawtooth wave in a firing range. Again, you can listen to it in chap10.ipynb.

图 10-5显示了此计算的示意图，其中f是锯齿波，g是脉冲响应，h是g的移位、缩放副本之和。

Figure 10-5 shows a diagram of this computation, where f is the sawtooth, g is the impulse response, and h is the sum of the shifted, scaled copies of g.

图 10-5. g 的缩放和平移副本之和的示意图。

例如：

For the example shown:

更广泛地说：

And more generally:

你可能在“卷积”中见过这个等式。它是f和g的卷积。这表明，如果输入是f，系统的脉冲响应是g ，则输出是f和g的卷积。

You might recognize this equation from “Convolution”. It is the convolution of f and g. This shows that if the input is f and the impulse response of the system is g, the output is the convolution of f and g.

总而言之，有两种方式可以思考系统对信号的影响：

In summary, there are two ways to think about the effect of a system on a signal:

1. 输入是一系列脉冲，因此输出是脉冲响应的缩放、移位副本之和；该和是输入和脉冲响应的卷积。

2. The input is a sequence of impulses, so the output is the sum of scaled, shifted copies of the impulse response; that sum is the convolution of the input and the impulse response.

3. 脉冲响应的离散傅里叶变换 (DFT) 是一个传递函数，它将系统对每个频率分量的影响编码为幅值和相位偏移。输入信号的 DFT 编码了其包含的频率分量的幅值和相位偏移。将输入信号的 DFT 与传递函数相乘即可得到输出信号的 DFT。

4. The DFT of the impulse response is a transfer function that encodes the effect of the system on each frequency component as a magnitude and phase offset. The DFT of the input encodes the magnitude and phase offset of the frequency components it contains. Multiplying the DFT of the input by the transfer function yields the DFT of the output.

这些描述的等价性并不令人意外。它本质上是对卷积定理的表述：时域中函数f和g的卷积对应于频域中的乘法。

The equivalence of these descriptions should not be a surprise. It is basically a statement of the Convolution Theorem: convolution of f and g in the time domain corresponds to multiplication in the frequency domain.

如果你想知道为什么卷积的定义是这样的（这在我们讨论平滑和差分窗时似乎有点反常），现在你知道原因了：卷积的定义自然地出现在 LTI 系统对信号的响应中。

And if you wondered why convolution is defined as it is, which seemed backward when we talked about smoothing and difference windows, now you know the reason: the definition of convolution appears naturally in the response of an LTI system to a signal.

卷积定理的证明

Proof of the Convolution Theorem

好了，我已经拖延够久了。是时候证明卷积定理（CT）了，该定理指出：

Well, I’ve put it off long enough. It’s time to prove the Convolution Theorem (CT), which states:

其中f和g是长度相同的向量，N。

where f and g are vectors with the same length, N.

我将分两步进行：

I’ll proceed in two steps:

1. 我将证明，在f为复指数的特殊情况下，与g进行卷积的效果是将f乘以一个标量。

2. I’ll show that in the special case where f is a complex exponential, convolution with g has the effect of multiplying f by a scalar.

3. 在更一般的情况下，如果f不是复指数函数，我们可以使用 DFT 将其表示为指数分量之和，计算每个分量的卷积（通过乘法），然后将结果相加。

4. In the more general case where f is not a complex exponential, we can use the DFT to express it as a sum of exponential components, compute the convolution of each component (by multiplication), and then add up the results.

这些步骤共同证明了卷积定理。但首先，让我们把需要的材料组装起来。函数g的离散傅里叶变换（我称之为G）为：

Together these steps prove the Convolution Theorem. But first, let’s assemble the pieces we’ll need. The DFT of g, which I’ll call G, is:

其中k为频率索引，取值范围为 0 到n ，n为时间索引，取值范围为 0 到n。结果是一个包含N 个复数的向量。

where k is an index of frequency from 0 to and n is an index of time from 0 to . The result is a vector of N complex numbers.

F的逆 DFT（我将其记为f）为：

The inverse DFT of F, which I’ll call f, is:

以下是卷积的定义：

Here’s the definition of convolution:

其中m是另一个时间索引，取值范围为 0 到。卷积运算满足交换律，因此我可以等价地写成：

where m is another index of time from 0 to . Convolution is commutative, so I could equivalently write:

现在让我们考虑f为频率为k 的复指数函数的特殊情况，我将其称为e _k：

Now let’s consider the special case where f is a complex exponential with frequency k, which I’ll call e_k:

其中k为频率索引，n为时间索引。

where k is an index of frequency and n is an index of time.

将e _k代入卷积的第二个定义，得到：

Plugging e_k into the second definition of convolution yields:

我们可以将第一个术语拆分成一个乘积：

We can split the first term into a product:

前半部分与m无关，因此我们可以将其从求和式中分离出来：

The first half does not depend on m, so we can pull it out of the summation:

现在我们知道第一项是e _k，求和式为（用m作为时间索引）。因此我们可以写成：

Now we recognize that the first term is e_k, and the summation is (using m as the index of time). So we can write:

这表明，对于每个复指数e _k，与g的卷积相当于将e _k乘以。用数学术语来说，每个e _k都是此运算的特征向量，而是相应的特征值（参见“微分”）。

which shows that for each complex exponential, e_k, convolution with g has the effect of multiplying e_k by . In mathematical terms, each e_k is an eigenvector of this operation, and is the corresponding eigenvalue (see “Differentiation”).

现在来看证明的第二部分。如果输入信号f不是复指数函数，我们可以通过计算其 DFT F将其表示为复指数函数的和。对于从 0 到 的每个k值，是频率为k的分量的复幅度。

Now for the second part of the proof. If the input signal, f, doesn’t happen to be a complex exponential, we can express it as a sum of complex exponentials by computing its DFT, F. For each value of k from 0 to , is the complex magnitude of the component with frequency k.

根据证明的第一部分，每个输入分量都是幅度为的复指数，因此每个输出分量都是幅度为的复指数。

Each input component is a complex exponential with magnitude , so each output component is a complex exponential with magnitude , based on the first part of the proof.

由于该系统是线性的，因此输出就是各个输出分量的总和：

Because the system is linear, the output is just the sum of the output components:

代入e _k的定义式，得到：

Plugging in the definition of e_k yields:

等式右侧是乘积的逆密度泛函理论(DFT) 。因此：

The righthand side is the inverse DFT of the product . Thus:

替换和：

Substituting and :

最后，对等式两边进行离散傅里叶变换，即可得到卷积定理：

Finally, taking the DFT of both sides yields the Convolution Theorem:

证毕。

QED.

练习

Exercises

这些练习的答案在chap10soln.ipynb……

Solutions to these exercises are in chap10soln.ipynb.

脉冲卷积

Convolution with Impulses

正如我们在“系统与卷积”中看到的那样，将一个信号与一系列脉冲进行卷积，其效果是将该信号的移位、缩放副本相加。

As we saw in “Systems and Convolution”, convolution of a signal with a series of impulses has the effect of adding up shifted, scaled copies of the signal.

举例来说，我会读取一个听起来像哔哔声的信号：

As an example, I’ll read a signal that sounds like a beep:

文件名 = '253887__themusicalnomad__positive-beeps.wav'
wave = thinkdsp.read_wave(filename)
wave.normalize()

filename = '253887__themusicalnomad__positive-beeps.wav'
wave = thinkdsp.read_wave(filename)
wave.normalize()

我将构建一个包含四个脉冲的波：

And I’ll construct a wave with four impulses:

imp_sig = thinkdsp.Impulses([0.005, 0.3, 0.6, 0.9],
                       amps=[1, 0.5, 0.25, 0.1])
impulses = imp_sig.make_wave(start=0, duration=1.0,
                             帧率=wave.帧率)

imp_sig = thinkdsp.Impulses([0.005, 0.3, 0.6,  0.9], 
                       amps=[1,     0.5, 0.25, 0.1])
impulses = imp_sig.make_wave(start=0, duration=1.0, 
                             framerate=wave.framerate)

然后将它们进行卷积：

and then convolve them:

卷积 = wave.convolve(脉冲)

convolved = wave.convolve(impulses)

图 11-1显示了结果，左上角为信号，左下角为脉冲，右侧为结果。

Figure 11-1 shows the results, with the signal in the top left, the impulses in the lower left, and the result on the right.

图 11-1.将一个信号（左上）与一系列脉冲（左下）进行卷积的效果。结果（右）是信号的平移和缩放副本之和。

你可以听听结果chap11.ipynb；听起来像是一连串音量逐渐减小的四声哔哔声。

You can hear the result in chap11.ipynb; it sounds like a series of four beeps with decreasing loudness.

这个例子的目的在于说明脉冲卷积会生成平移缩放后的副本。这个结果在下一节中会很有用。

The point of this example is just to demonstrate that convolution with impulses makes shifted, scaled copies. This result will be useful in the next section.

调幅

Amplitude Modulation

幅度调制（AM）技术用于广播调幅广播等多种用途。在发射机端，信号（可能包含语音、音乐等）通过与余弦信号（作为“载波”）相乘进行“调制”。调制后的信号为适合无线电广播的高频波。美国调幅广播的典型频率范围为 500–1600 kHz（参见https://en.wikipedia.org/wiki/AM_broadcasting）。

Amplitude modulation (AM) is used to broadcast AM radio, among other applications. At the transmitter, the signal (which might contain speech, music, etc.) is “modulated” by multiplying it with a cosine signal that acts as a “carrier wave”. The result is a high-frequency wave that is suitable for broadcast by radio. Typical frequencies for AM radio in the United States are 500–1600 kHz (see https://en.wikipedia.org/wiki/AM_broadcasting).

在接收端，广播信号被“解调”以恢复原始信号。令人惊讶的是，解调的工作原理是将广播信号再次乘以同一个载波。

At the receiving end, the broadcast signal is “demodulated” to recover the original signal. Surprisingly, demodulation works by multiplying the broadcast signal, again, by the same carrier wave.

为了了解其工作原理，我将用一个频率为 10 kHz 的载波来调制一个信号。信号如下：

To see how that works, I’ll modulate a signal with a carrier wave at 10 kHz. Here’s the signal:

文件名 = '105977__wcfl10__favorite-station.wav'
wave = thinkdsp.read_wave(filename)
wave.unbias()
wave.normalize()

filename = '105977__wcfl10__favorite-station.wav'
wave = thinkdsp.read_wave(filename)
wave.unbias()
wave.normalize()

这就是承运商：

And here’s the carrier:

rier_sig = thinkdsp.CosSignal(freq=10000)
carrier_wave = carrier_sig.make_wave(duration=wave.duration,
                                     帧率=wave.帧率)

carrier_sig = thinkdsp.CosSignal(freq=10000)
carrier_wave = carrier_sig.make_wave(duration=wave.duration, 
                                     framerate=wave.framerate)

我们可以使用运算符将​​它们相乘，该运算符逐元素地*乘以波数组：

We can multiply them using the * operator, which multiplies the wave arrays elementwise:

调制波 = 波 * 载波

modulated = wave * carrier_wave

结果听起来很糟糕。你可以听听chap11.ipynb。

The result sounds pretty bad. You can hear it in chap11.ipynb.

图 11-2显示了频域中的情况。第一行是原始信号的频谱。第二行是乘以载波后的调制信号的频谱。它包含原始频谱的两个副本，分别偏移了正负 10 kHz。

Figure 11-2 shows what’s happening in the frequency domain. The top row is the spectrum of the original signal. The next row is the spectrum of the modulated signal, after multiplying by the carrier. It contains two copies of the original spectrum, shifted by plus and minus 10 kHz.

为了理解这一点，需要记住，时域中的卷积对应于频域中的乘法。反之亦然，时域中的乘法对应于频域中的卷积。当我们将信号乘以载波时，实际上是将信号的频谱与载波的离散傅里叶变换（DFT）进行卷积。

To understand why, recall that convolution in the time domain corresponds to multiplication in the frequency domain. Conversely, multiplication in the time domain corresponds to convolution in the frequency domain. When we multiply the signal by the carrier, we are convolving its spectrum with the DFT of the carrier.

由于载波是一个简单的余弦波，其离散傅里叶变换 (DFT) 为两个脉冲，频率分别为±10 kHz。与这两个脉冲进行卷积运算，即可得到频谱的移位和缩放副本。请注意，调制后频谱的幅度变小了。原始信号的能量被分配到这两个副本中。

Since the carrier is a simple cosine wave, its DFT is two impulses, at plus and minus 10 kHz. Convolution with these impulses makes shifted, scaled copies of the spectrum. Notice that the amplitude of the spectrum is smaller after modulation. The energy from the original signal is split between the copies.

我们再次乘以载波来解调信号：

We demodulate the signal by multiplying by the carrier wave again:

解调后的值 = 调制后的值 * 载波

demodulated = modulated * carrier_wave

图 11-2的第三行显示了结果。同样，时域中的乘法对应于频域中的卷积，从而生成频谱的平移和缩放副本。

The third row of Figure 11-2 shows the result. Again, multiplication in the time domain corresponds to convolution in the frequency domain, which makes shifted, scaled copies of the spectrum.

图 11-2.幅度调制演示。第一行是信号的频谱；第二行是调制后的频谱；第三行是解调后的频谱；最后一行是经过低通滤波后的解调信号。

由于调制后的频谱包含两个峰值，每个峰值都被分成两半，并分别向正负 20 kHz 方向移动。其中两份副本在 0 kHz 处相遇并叠加；另外两份副本最终分别以正负 20 kHz 为中心。

Since the modulated spectrum contains two peaks, each peak gets split in half and shifted by plus and minus 20 kHz. Two of the copies meet at 0 kHz and get added together; the other two copies end up centered at plus and minus 20 kHz.

如果你听一下解调后的信号，会发现它听起来相当不错。额外的频谱副本增加了原始信号中没有的高频成分。这些高频成分非常高，大多数扬声器无法播放，大多数人也听不到，但如果你拥有好的扬声器和灵敏的耳朵，或许就能听到。

If you listen to the demodulated signal, it sounds pretty good. The extra copies of the spectrum add high-frequency components that were not in the original signal. These are so high that most speakers can’t play them and most people can’t hear them, but if you have good speakers and good ears, you might.

在这种情况下，您可以通过应用低通滤波器来去除多余的成分：

In that case, you can get rid of the extra components by applying a low-pass filter:

解调频谱 = 解调.生成频谱(full=True)
解调频谱.低通(10000)
过滤后的频谱 = 解调频谱.生成波形()

demodulated_spectrum = demodulated.make_spectrum(full=True)
demodulated_spectrum.low_pass(10000)
filtered = demodulated_spectrum.make_wave()

结果与原始波形非常接近，尽管解调和滤波后损失了大约一半的功率。这在实际应用中不是问题，因为广播信号的发射和接收过程中会损失更多的功率。既然无论如何我们都需要放大结果，那么再放大两倍也无关紧要。

The result is quite close to the original wave, although about half of the power is lost after demodulating and filtering. That’s not a problem in practice, because much more of the power is lost in transmitting and receiving the broadcast signal. Since we have to amplify the result anyway, another factor of 2 is not an issue.

采样

Sampling

我解释幅度调制，一部分原因是它很有趣，但更主要的原因是它有助于我们理解采样。采样是指在一系列时间点（通常间隔相等）测量模拟信号的过程。

I explained amplitude modulation in part because it is interesting, but mostly because it will help us understand sampling. Sampling is the process of measuring an analog signal at a series of points in time, usually with equal spacing.

例如，我们用作示例的 WAV 文件是通过使用模数转换器 (ADC) 对麦克风输出进行采样而录制的。它们大多采用 44.1 kHz 的采样率，这是“CD 音质”的标准采样率；或者采用 48 kHz 的采样率，这是 DVD 音质的标准采样率。

For example, the WAV files we have used as examples were recorded by sampling the output of a microphone using an analog-to-digital converter (ADC). The sampling rate for most of them is 44.1 kHz, which is the standard rate for “CD-quality” sound, or 48 kHz, which is the standard for DVD sound.

在 48 kHz 时，折叠频率为 24 kHz，这比大多数人能听到的频率要高（参见https://en.wikipedia.org/wiki/Hearing_range）。

At 48 kHz, the folding frequency is 24 kHz, which is higher than most people can hear (see https://en.wikipedia.org/wiki/Hearing_range).

在这些波形中，每个样本有 16 位，因此有 2^ ¹⁶个不同的级别。这种“位深度”已经足够，增加更多位数并不会明显改善音质（参见https://en.wikipedia.org/wiki/Digital_audio）。

In most of these waves, each sample has 16 bits, so there are 2¹⁶ distinct levels. This “bit depth” turns out to be enough that adding more bits does not improve the sound quality noticeably (see https://en.wikipedia.org/wiki/Digital_audio).

当然，除了音频信号之外的其他应用可能需要更高的采样率来捕捉更高的频率，或者更高的位深度来更保真地再现波形。

Of course, applications other than audio signals might require higher sampling rates in order to capture higher frequencies, or higher bit depth in order to reproduce waveforms with more fidelity.

为了演示采样过程的效果，我将从一个以 44.1 kHz 采样率采集的波形开始，并从中选取大约 11 kHz 的样本。这与从模拟信号中采样并不完全相同，但效果是一样的。

To demonstrate the effect of the sampling process, I am going to start with a wave sampled at 44.1 kHz and select samples from it at about 11 kHz. This is not exactly the same as sampling from an analog signal, but the effect is the same.

首先，我将加载一段鼓独奏的录音：

First I’ll load a recording of a drum solo:

文件名 = '263868__kevcio__amen-break-a-160-bpm.wav'
wave = thinkdsp.read_wave(filename)
wave.normalize()

filename = '263868__kevcio__amen-break-a-160-bpm.wav'
wave = thinkdsp.read_wave(filename)
wave.normalize()

图 11-3（上）显示了该波的频谱。下面是对该波进行采样的函数：

Figure 11-3 (top) shows the spectrum of this wave. Now here’s the function that samples from the wave:

def sample(wave, factor=4):
    ys = np.zeros(len(wave))
    ys[::factor] = wave.ys[::factor]
    返回 thinkdsp.Wave(ys, framerate=wave.framerate)

def sample(wave, factor=4):
    ys = np.zeros(len(wave))
    ys[::factor] = wave.ys[::factor]
    return thinkdsp.Wave(ys, framerate=wave.framerate)

我将用它来选择每四个元素中的一个：

I’ll use it to select every fourth element:

sampled = sample(wave, 4)

sampled = sample(wave, 4)

结果帧速率与原信号相同，但大部分元素的值都为零。如果播放采样后的波形，听起来不太好。采样过程引入了原信号中原本不存在的高频成分。

The result has the same frame rate as the original, but most of the elements are zero. If you play the sampled wave, it doesn’t sound very good. The sampling process introduces high-frequency components that were not in the original.

图 11-3（底部）显示了采样波的频谱。它包含原始频谱的四个副本（看起来像五个副本，因为其中一个副本被分割在最高频率和最低频率之间）。